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not  include  condensation  processes ,  but  it  does  include  a  radiation  pnrareterizat; 
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moisture  between  the  boundnrv  Inver  and  the  earth's  sur'acc.  i  mil  at  ions  have  bet 
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1 .  Introduction 


Trajectory  and  advection  models  (Muench  1983;  Muench  and  Chisholm  1985; 
Muench  1989)  and  automated  display  systems  such  as  McIDAS  and  PROFS  (Schlatter 
et  al.  1985)  show  promise  of  providing  added  guidance  for  the  very  short  term 
forecasting  range  of  0  to  6  hours.  Full  synoptic  models  such  as  the  LFM  and  NGM 
have  skill  from  about  12  hours  to  beyond  48  hours.  However,  both  of  these 
approaches  fall  short  in  the  6  to  18  hour  period  that  is  perhaps  the  most  vital  to 
the  terminal  forecaster.  It  is  in  this  period  that  meso-/9  scale  disturbances  (such  as 
fronts  or  squall  lines)  that  are  in  the  local  area  will  have  an  immediate  impact  on 
the  terminal  forecast.  Advection  models  begin  failing  after  a  few  hours  because  of 
synoptic  and  mesoscale  changes  in  wind  patterns  and  because  topographic  and 
geographic  influences  are  often  not  incorporated  (Muench  1983;  Muench  1989).  On 
the  other  hand,  the  initialization  of  synoptic  scale  models  is  based  on  synoptic  scale 
data  that  does  not  normally  include  mesoscale  features  present  at  the  time  of 
initialization,  and  the  coarse  grid  spacing  does  not  allow  mesoscale  features  to 
develop  during  the  integration.  Hence,  the  synoptic  scale  model  cannot  be  expected 
to  forecast  these  mesoscale  features.  Other  studies,  however,  have  shown  thai. 
mesoscale  models  can  forecast  the  mesoscale  structure  of  a  variety  of  atmospheric 
phenomena  even  when  initialized  with  observational  data  of  only  somewhat  finer 
resolution  than  routinely  used  in  synoptic  scale  models  (Zhang  and  Fritsch  1986;  and 
others). 

There  are  currently  many  mesoscale  numerical  models  within  the  research 
community  Isee  Keyser  and  Uccellini  (1987)  for  a  discussion  of  some  of  these 
models].  Each  of  these  models  was  developed  initially  to  look  at  certain  subsets  of 
meteorological  phenomena  and  then  evolved  into  much  more  complex  modeling 
systems  capable  of  simulating  a  wide  range  of  atmospheric  conditions.  While  these 
models  differ  substantially  in  details,  such  as  the  level  of  complexity  for  the 
microphysical  parameterizations  or  turbulence  closure,  they  share  some  basic 
features.  Most  are  what  can  be  considered  “research  models”  in  that  they  are  not 


designed  specifically  for  operational  use;  their  primary  use  is  to  try  to  evaluate  the 
role  of  various  atmospheric  processes  for  better  understanding  of  the  dynamics  and 
physics  of  the  atmosphere  on  the  mesoscale.  In  general,  these  models  tend  to 
include  the  most  realistic  boundary  layer,  radiation,  and  microphysical 
parameterizations  possible  and  still  achieve  acceptable  computation  times  on  the 
current  generation  of  supercomputers.  They  are  usually  composed  of  at  least  one 
level  of  horizontal  nesting  and  employ  enough  vertical  layers  to  resolve  explicitly 
the  evolving  boundary  layer.  Recently,  the  ETA  mesoscale  model  (Mesinger  et  al. 
1988;  Mesinger  et  al.  1990)  has  been  used  in  an  operational  mode  (Kalnay  et  al.  1991), 
and  a  simplified  version  of  the  PSL'/NCAR  model  has  been  used  in  a  pseudo- 
operational  setting  (Warner  and  Seaman  1990).  These  models  are  somewhat  simpler 
than  the  current  generation  of  research  models,  but  the  ETA  model  still  requires  a 
supercomputer,  and  the  simplified  PSU/NCAR  model  requires  significant 
computational  time  on  a  conventional  mainframe. 

The  model  described  here  represents  a  significant  deviation  from  this 
approach.  Here  we  attempt  to  develop  a  model  specifically  for  operational  use  on 
the  current  or  near-future  generation  of  super-micro  computers  by  making  use  of  a 
consistent  set  of  approximations  and  parameterizations  that  incorporates  the 
physical  understanding  of  mesoscale  processes  emerging  from  the  results  of  research 
models.  Our  goal  is  a  mesoscale  model  that  is  much  faster  computationally  than 
research  models  but  still  has  enough  skill  to  provide  guidance  to  the  terminal 
forecaster  in  the  6—18  h  forecast  window.  We  seek  to  develop  a  mesoscale  model 
that  can  be  run  locally  by  each  forecast  office,  centered  on  its  location,  to  provide 
local  guidance. 

Given  the  rapid  increase  in  speed  and  power  of  small  computers,  it  could  be 
argued  that  the  current  research  models  could  be  run  locally  on  relatively 
inexpensive  computers  within  the  next  decade.  There  is  currently  some  debate 
within  the  meteorological  community  on  the  relative  merits  of  running  highly 
complex  mesoscale  models  locally  at  each  forecast  office  compared  to  a  centralized 
guidance  with  distributed  products  similar  to  those  provided  for  synoptic  scale 
models  (McPherson  1991).  Many  argue  that  a  very  complex  mesoscale  model 
requiring  a  major  data  assimilation  effort  and  excellent  boundary  conditions 
produced  by  synoptic  scale  models  can  only  be  maintained  effectively  at  a  national 
center.  Further,  these  individuals  argue  that  running  a  near-research-grade  model 


at  many  sites  with  largely  overlapping  domains  is  a  waste  of  computer  resources. 
This  work  supports  this  philosophy  of  mesoscale  guidance,  but  argues  that 
simplified,  computationally  fast,  mesoscale  models  can  provide  a  valued-added 
service  when  run  at  local  forecast  offices.  While  the  simple  model  described  in  this 
report  cannot  be  expected  to  give  as  skillful  guidance  over  a  forecast  period  as  a 
research  grade  model,  it  can  be  run  locally  long  before  the  more  complex  model 
output  is  available,  providing  a  “first  look” — similar  to  the  way  the  National 
Meteorological  Center’s  LFM  model  is  used  operationally  today  (Petersen  and 
Stackpole  1989).  With  future  advances,  a  model  of  the  type  developed  here  could 
possibly  be  run  in  an  interactive  fashion-similar  to  interactive  sounding 
modification  programs  now  available  that  allow  the  forecaster  to  investigate  “what 
if”  scenarios. 

While  it  is  possible  to  produce  a  simplified  version  of  a  research  model  that 
is  computationally  efficient  enough  to  provide  operational  mesoscale  forecasts 
(Warner  and  Seaman  1990),  we  have  taken  a  different  approach.  Several  studies 
have  identified  topographic  forcing  and  boundary-layer  processes  as  major 
contributors  to  mesoscale  structure  (Gauntlet  et  al.  1984;  Ookouchi  et  al.  1984; 
Beniamin  and  Carlson  1986;  Nickerson  et  al.  1986;  and  others),  so  any  successful 
mesoscale  model  must  include  these  processes.  The  sigma-coordinate  system  used  in 
many  models  allows  topography  to  be  included  in  a  straightforward  manner.  Most 
models  use  many  layers  near  the  surface  to  resolve  boundary-layer  processes,  but 
this  adds  to  the  computational  expense  by  increasing  the  number  of  grid  points  in 
the  model.  Our  approach  is  to  treat  the  boundary  layer  as  a  single  layer  of  known 
structure  that  can  evolve  over  the  course  of  the  integration.  The  model  equations 
are  recast  in  a  coordinate  system  based  on  the  depth  of  the  evolving  boundary 
layer.  This  allows  the  boundary-layer  structure  and  fluxes  to  be  represented  in  a 
way  similar  to  that  used  in  some  mixed-layer  models  (Lavoie  1972;  Colby  1984;  and 
others)  so  that  the  physical  processes  are  represented  without  the  computational 
expense  associated  with  a  large  number  of  vertical  layers. 

This  paper  describes  a  prototype  dry  mesoscale  model  that  has  been 
developed  using  this  approach.  Section  2  introduces  the  boundary-layer  coordinate 
system  on  which  the  mode!  is  based  and  presents  the  numerical  details  of  the  model. 
Section  3  summarizes  the  boundary  layer  and  radiation  parameterizations  employed 
in  the  model.  This  is  followed,  in  section  4,  by  the  results  of  simulations  designed 
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to  assess  the  capabilities  of  the  model.  Section  5,  gives  the  results  of  some  very 
preliminary  tests  using  complex  terrain  and  describes  how  the  terrain  field  must  be 
adjusted  to  be  consistent  with  the  boundary  conditions.  The  paper  concludes  with 
discussion  and  plans  for  extensions  to  the  dry  model. 


2.  Model  description 

a.  Basic  model  equations 

As  discussed  in  the  Introduction,  most  mesoscale  models  rely  on  a  large 
number  of  layers  near  the  surface  in  order  to  resolve  explicitly  the  growth  and 
decay  of  the  planetary  boundary  layer  and  to  simulate  correctly  the  fluxes  of  heal 
and  moisture  that  couple  the  atmosphere  to  the  surface.  This  leads  to  models  with 
a  large  number  of  layers  in  the  vertical,  resulting  in  large  numbers  of  computations 
per  timestep  and  relatively  long  computation  times  even  in  supercomputer 
environments  (Peilke  1984;  his  Appendix  B).  On  the  other  hand,  some  early  models 
showed  success  in  simulating  nonboundary-layer-driven  flows  (such  as  mountain  lee- 
waves)  with  relatively  few  layers  (Anthes  and  Warner  1978),  and  other  models 
looking  specifically  at  boundary-layer  processes  were  able  to  capture  successfully 
the  important  features  by  treating  the  boundary  layer  as  a  single  layer  which  could 
dynamically  grow  and  collapse  (Lavoie  1972;  Keyser  and  Anthes  1977;  Anthes  et  al. 
1982;  Colby  1984). 

The  present  formulation  seeks  to  marry  these  two  approaches  into  a  single 
three-dimensional  mesoscale  model.  The  lowest  layer  of  the  model  is  the  boundary 
layer-  it  is  allowed  to  grow  in  depth  or  collapse  at  each  grid  point  as  the  simulation 
proceeds  in  response  to  changes  produced  by  the  boundary-layer  and  radiation 
parametenzations.  The  layers  above  the  boundary  layer  adjust  dynamically  to  the 
changing  boundary-layer  depth,  while  simulating  the  horizontal  and  vertical 
advections  of  momentum  and  thermodynamic  variables  and  maintaining  the  proper 
balances  that  hold  at  the  meso-/3  scale. 

The  vertical  coordinate  that  allows  for  this  changing  boundary-layer  depth  is 
a  modification  of  the  common  o-coordinate  [cr  =(p  —  pt) /(ps  —  pt )  where  ps  is  the 
surface  pressure  and  pt  is  a  pressure  level  specified  as  the  top  of  the  model  (we 
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take  pt  =  100  mb  for  this  study)!.  If  we  let  ah  represent  the  height  of  the  top  of 
the  boundary  layer  in  a  (see  Fig.l),  we  can  define  a  new  vertical  coordinate,  f],  as 


T] 


H  = 


On 

1  —  <T„ 


o  <  ah 
o  >  ah 


(2.1) 


We  refer  to  the  V  coordinate  system  as  “boundary-layer  coordinates”  [note  that  this 
is  not  related  to  what  have  been  called  ETA  models  (Mesinger  et  al.  1988;  Mesmger 
et  al.  1990)]  The  vector  momentum  equation,  hydrostatic  relation,  continuity 
equation,  thermodynamic  equation,  and  specific  humidity  conservation  equation  can 
be  written  in  the  T]  system  as  (see  appendix  B  for  details  on  the  derivation  of  this 
set) 
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where  ir  =  ps—pt,  and  the  other  terms  have  their  normal  meteorological  definitions 
(see  Appendix  A  for  a  list  of  variables).  It  should  be  noted  that  the  a  which 
appears  on  the  right-hand  side  of  (2.2)  is  a  dependent  variable  which,  in  general,  will 
not  be  constant  on  a  constant-^  surface.  Expansion  of  this  gradient  of  a  product 
results  in  an  additional  pressure-gi  adient-force  type  term  not  present  in  traditional 
cr-coordinate  models.  It  should  also  be  noted  that  the  humidity  is  treated  here  as  a 
passive  scalar,  so  (2.6)  simply  represents  conservation  of  q  (that  is,  dq/dt  =  0)  in 
the  r/  system.  When  condensation  and  evaporation  are  included,  appropriate  source 
and  sink  terms  need  to  be  added  to  the  rhs  of  (2.6)  and  another  conservation 
equation  governing  the  condensate  is  required. 
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Fig.  1.  Schematic  structure  of  the  five-layer  ^-coordinate  model. 


Equation  (2.4)  is  not  used  directly  in  the  model,  but  it  does  provide  the 
means  of  calculating  vertical  velocities  and  the  rate  of  change  of  surface  pressure. 
Integration  of  (2.4)  from  the  model  top  to  the  surface  yields 

If  -  -  f  \&H 

(2.7) 

-r  ~iH -rv)\dv 

a:r  dy  j 

where  H-  —  1  —  cr^  and  H.  —  cr>,  represent  the  values  of  H  above  and  below  the 
boundary-layer  top,  respectively,  as  given  by  (2.1).  Integration  of  (2.4)  from  the 
top  of  the  model  down  to  a  specific  interface  level,  along  with  the  use  of  dr  /dt 
found  with  (2.7)  allows  the  determination  of  (Hi 7)  at  each  interface  as 
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(Hr?)  =  -(77  +  1) 


dr?  .  (2.8) 


+  "‘I]  -  Zi&H'ru)  +  A(H‘™ 

It  is  not  difficult  to  show  that  the  continuous  equations  in  boundary-layer 
coordinates  satisfy  the  proper  energy  conservation  constraints  (Colby  and  Seitter, 
1990).  It  is  desirable  to  have  the  finite  difference  forms  of  the  equations  satisfy 
these  energy  conservation  constraints  as  well,  which  can  provide  guidance  on  the 
form  the  finite  difference  equations  should  take.  In  particular,  the  energy 
constraints  lead  to  specific  forms  for  the  solution  of  the  thermodynamic  equation, 
(2.5),  and  the  integration  of  the  hydrostatic  equation,  (2.3),  to  find  the  geopotentials 
of  midlevels  of  the  Tj- layers.  The  resulting  finite  difference  forms  are  considerably 
different  in  structure  from  equations  (2.2)—(2.6),  and  are  given  in  Appendix  B. 

The  last  terms  on  the  rhs  of  (2.2),  and  (2.5)  include  a  “friction”  term,  F, 
which,  above  the  boundary  layer,  is  given  by  a  horizontal  eddy  diffusion.  This 
term  is  modeled  by  a  simple  Fickian  diffusion,  KV2#,  where  K  is  a  constant  eddy 
viscosity  and  ^  is  the  variable  of  interest  hi,  v,  or  T). 

b.  Grid  domain  and  horizontal  nesting 

A  staggered  grid  is  used  in  both  the  vertical  and  horizontal  directions.  In 
the  vertical,  all  variables  are  layer  quantities  except  vertical  velocities,  which  are 
defined  at  interface  levels  (see  Fig.  1).  Although  the  coordinate  system  is  quite 
general  and  would  allow  multiple  layers  in  the  boundary  layer  (1  >  r?  >  0),  we 
currently  run  the  model  with  the  boundary  layer  representing  a  single  model  layer 
and  specify  four  layers  above  the  boundary  layer  that  are  equally  spaced  in  V. 
Horizontally,  velocities  are  defined  on  staggered  points  that  surround  the  points  on 
which  all  other  variables  are  defined,  as  in  Anthes  and  Warner  (1978).  In  order  to 
increase  the  overall  model  domain  size  and  move  the  lateral  boundaries  away  from 
the  area  of  primary  interest,  a  horizontal  nesting  of  the  model  is  employed  as 
developed  by  Zhang  et  al.  (1986).  A  fine  grid  mesh  (FGM)  with  20  km  resolution  is 
nested  in  a  coarse  grid  mesh  (CGM)  with  60  km  resolution.  The  complete  nested 
staggered  grid  domain,  showing  all  points,  is  shown  in  Fig.  2.  A  3:1  ratio  of  FGM 
points  to  CGM  points  is  necessary  with  a  staggered  grid  so  that  both  “velocity”  (x’s 
in  Fig.  2)  and  “thermodynamic”  (o's  in  Fig.  2)  points  can  be  coincident  in  the  overlap 
region  (Zhang  et  al.  1986).  The  CGM  domain  covers  1320  km  x  1320  km  while  the 
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Fig.  2.  Horizontal  plot  of  the  nested  staggered  grid  structure 
used  in  the  model.  Velocities  are  specified  on  all  ‘V  points 
and  all  other  variables  are  defiend  on  “o”  points  (referred  to 
as  “thermodynamic”  points  in  the  text. 
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FGM  domain  is  480  km  x  480  km  for  “thermodynamic”  points  (all  displays  will  be 
made  on  “thermodynamic”  point  arrays,  with  any  displayed  velocities  being  averaged 
to  these  points). 

The  two-way  interactive  nesting  procedure  of  Zhang  et  J.  (1986)  is  used  with 
a  few  minor  modifications.  In  the  calculation  of  tendencies  for  the  “thermodynamic" 
points  in  the  FGM,  for  example,  a  simple  linear  interpolation  is  used  between  CGM 
points  nearest  the  boundary  FGM  point  rather  than  the  “Lagrangian  interpolation” 
used  by  Zhang  et  al.  (1986).  Also,  rather  than  use  a  Newtonian  damping  scheme  near 
the  boundaries,  the  eddy  diffusion  coefficient  K  is  increased  for  FGM  points  within 
two  grid  intervals  of  the  FGM/CGM  interface.  This  increased  diffusion  is  applied 
only  in  the  momentum  equation  to  help  control  noise  resulting  from  the 
overspecification  there,  and  has  a  net  effect  quite  similar  to  the  Newtonian  damping 
term  used  by  Zhang  et  al.  (see  Colby  and  Seitter  1990;  Kurihara  and  Bender  1980). 

c.  Time  integration  and  outer  lateral  boundary  condition 

Time  integration  for  the  model  is  performed  using  the  leapfrog  scheme  with 
an  Asselin  filter  except  for  the  moisture  equation  (2.6),  which  employs  a  forward 
timestepping  scheme  with  upstream  differencing.  The  Asselin  filter  factor  is  set  to 
0.5  as  suggested  by  Schlesinger  et  al.  (1983)  for  models  employing  moderate 
diffusion.  The  time  step  for  FGM  points  is  20  s  and  for  CGM  points  is  60  s.  The 
flow  relaxation  condition  of  Davies  (1976)  is  used  on  the  lateral  boundaries  of  the 
CGM  following  the  work  of  Seitter  (1987)  who  found  that  this  condition  was  well- 
behaved,  provided  a  simple  means  of  allowing  external  information  to  be  introduced 
into  the  model,  and  did  not  require  the  smoothing  operator  necessary  in  the  Ferkey 
and  Kreitzberg  (1976)  sponge.  The  flow  relaxation  condition  is  normally  applied 
over  a  5  gridpoint  wide  region  near  the  boundary,  and  solutions  in  this  “relaxation 
region"  should  be  considered  modified.  The  flow  relaxation  condition  is  only  applied 
at  the  lateral  boundaries  of  the  CGM  domain,  however,  so  FGM  points  are  not 
significantly  influenced  by  the  relaxation  region. 

Experimentation  has  shown  that  application  of  the  condition  over  a  6 
gridpoint  wide  region  yields  a  slightly  improved  transition  between  the  interior  and 
boundary  areas,  even  though  the  weight  on  the  sixth  gridpoint  from  the  boundary  is 
very  small.  The  weights  on  the  staggered  gridpoints  are  set  so  that  they  smoothly 
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match  the  logarithmic  decrease  of  the  weighting  function  from  the  boundary  to  the 
interior. 


d.  Coupling  of  the  model  with  the  boundary -layer  parameterization 

The  basic  equations  that  make  up  the  boundary-layer  parameterization  are 
described  in  detail  in  section  3.  The  boundary-layer  package  provides  the  model 
with  the  fluxes  of  heat,  moisture,  and  momentum  from  the  surface,  and  calculates 
the  rate  of  change  of  boundary-layer  height,  which  is  fundamental  to  the  boundary- 
layer  coordinate  formulation.  It  also  provides  tendencies  for  the  ground  surface 
temperature  and  ground  wetness,  which  are  then  timestepped  using  a  leapfrog  scheme 
with  Asselin  filter  just  as  the  other  prognostic  variables  in  the  model.  A  summary 
of  these  quantities  is  shown  in  Table  1. 

It  is  important  to  note  that  the  boundary-layer  routine  diagnoses  the  current 
boundary-layer  depth  (hydrostatically)  from  the  current  cr*,  and  returns  a  time  rate 
of  change  of  boundary-layer  height,  3h/3f,  where  h  is  in  geometric  heigiit  above 
the  ground.  The  model  requires  both  the  rate  of  change  of  boundary-layer  height 
and  the  new  boundary-layer  height  in  terms  of  cr,  so  a  conversion  must  be  made. 

The  hydrostatic  relation  may  be  written  in  o-coordinates  as 


3cr 


— ra 


(2.9) 


which  can  be  integrated  from  the  surface  (a  =  1)  to  the  top  of  the  boundary  layer 
(ct  =cr*)  to  obtain 

r  cr* 

0*  --  <t>s  =  —  r  a  do  (2.10) 

where  is  the  geopotential  of  the  surface  and  <ph  is  the  geopotentiai  of  the 
boundary-layer  top.  With  the  excellent  approximation  that  (<t>h  -  <ps)  =  gh,  and 
making  use  of  the  mean  value  theorem,  we  can  obtain 


h 


cr*; 


(2.11) 


Taking  the  partial  derivative  with  respect  to  time  of  this  expression  and 


—  10  — 


Table  1.  Quantities  provided  to  the  model  by  the  boundary-layer 
parameterization. 


Variable  returned 

Description 

dh 

dt 

Rate  of  change  in  height  of  boundary  layer  top 

dTG 

dt 

Rate  of  change  of  ground  temperature 

dG  W 
dt 

Rate  of  change  of  ground  wetness 

Tkb 

Diagnosed  “surface”  temperature 

Qk.b 

Diagnosed  “surface”  mixing  ratio 

SH 

Sensible  heat  flux 

LH 

Latent  heat  flux 

GS 

Surface  soil  heat  flux 

NR 

Net  radiation  (incoming  shortwave  minus 
outgoing  IR) 

T  xy 

Surface  stress  (friction) 

rearranging  yields 

9^  _ QJdh  ,  (1  -  oh)dr  ,  (1  ~  Qh) 3a  n 

dt  radt  ^  '  *  dt  a  dt  u  ; 

Equation  (2.12)  provides  a  means  of  computing  the  rate  of  change  of  from  dh  dt 
provided  we  can  estimate  all  the  other  terms  on  the  rhs.  Experimentation  has 
shown  that  the  last  term  is  always  at  least  two  orders  of  magnitude  smaller  than 
the  first  two  terms  on  the  rhs,  so  that  term  is  dropped.  The  remaining  information 
is  readily  available  from  the  mode!  since  the  boundary-layer  temperature  can  be 
used  to  find  a  and  3x /3 1  can  be  computed  using  (2.7).  The  tendency  produced  by 
(2.12)  is  used  on  the  current  timestep  where  required  [for  example  in  (2.8)),  and  it  is 
also  used  to  calculate  the  new  value  of  ah  in  a  leapfrog  timestepping  scheme  (with 
Asselin  filter)  identical  to  that  used  for  the  other  prognostic  variables.  A 
horizontal  diffusion  of  the  form  KV'oh  is  applied  to  the  tendency  of  oh,  with  K  = 
2  X  10s  m2s~'.  This  weakly  couples  the  boundary-layer  height  to  the  surrounding 
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grid  points  and  suppresses  exceedingly  sharp  gradients  or  noise  in  the  boundary- 
layer  height  field. 

e.  Horizontal  diffusion 

After  experimentation  with  simulations  of  a  variety  of  atmospheric 
phenomena  (Colby  and  Seitter  1990),  the  eddy  diffusion  coefficient  was  set  at  K  =  2 
x  10s  m2s-1  for  the  momentum  components  in  the  FGM.  This  yields  a 
nondimensional  eddy  viscosity  of  0.01  which  is  somewhat  smaller  than  the  value 
used  by  Anthes  and  Warner  (1978)  in  their  six-layer  model.  The  eddy  diffusion 
coefficient  in  the  CGM  is  set  to  K  =  6  x  10s  m2s_1  for  the  momentum  components 
in  order  to  yield  the  same  nondimensional  eddy  viscosity  there.  The  diffusion  of 
temperature  uses  an  eddy  diffusion  coefficient  one  order  of  magnitude  smaller  than 
that  used  for  the  momentum  components. 


3.  Boundary-layer  parameterization 

The  boundary-layer  parameterization  is  composed  of  several  parts:  radiation; 
surface  energy  and  moisture  balance;  surface  drag;  surface  temperature  and  specific 
humidity;  and  boundary-layer  height  change.  Each  of  the  formulations,  described  in 
the  following  subsections,  is  designed  to  balance  computational  speed  with  physical 
accuracy.  In  most  cases,  the  physical  accuracy  is  restricted  by  the  limited  vertical 
resolution  of  the  model  rather  than  by  the  mathematical  formulation.  The 
boundary-layer  parameterization  package  has  been  tested  in  a  one-dimensional  form 
against  data  taken  during  the  O’Neill  boundary  layer  experiment  (Lettau  and 
Davidson  1957)  and  found  to  reproduce  the  evolution  of  the  boundary-layer 
structure  and  surface  temperature  and  moisture  quite  well  (Colby  and  Seitter  1990V 

a.  Radiation 

The  radiation  parameterization  is  taken  from  Katayama  (1972)  and  is  a  routine 
originally  designed  for  use  in  the  UCLA  general  circulation  model.  Although  the 
results  presented  in  this  paper  are  for  a  prototype  model  with  no  condensation 
processes,  the  radiation  parameterization  allows  for  cloud  layers,  and  the  following 
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discussion  describes  this  more  complete  form.  The  incident  radiation  and  infra-red 
(IR)  emission  are  calculated  separately.  The  model  incorporates  an  exponential  fit 
to  the  data  for  specific  humidity  to  allow  simple  integration  of  water  content.  The 
concentration  of  C02  is  included  in  a  fixed  form  based  on  experimental  data,  and  its 
contribution  is  then  a  constant. 

The  influx  of  radiation  is  computed  by  starting  with  the  solar  constant  and 
modifying  it  for  albedo  at  the  top  of  the  atmosphere.  Scattered  and  absorbable 
radiation  are  computed  separately,  the  fraction  being  assumed  constant  (35% 
available  for  absorption,  65%  scattered  to  the  ground).  The  scattered  part  of  the 
incident  radiation  is  corrected  for  multiple  reflection  between  the  atmosphere  and 
the  ground  and  given  by 


GLW*  =  (0.651  )S0cos(ZT)(j^-^-) 


(3.1) 


If  a  cloud  layer  is  present  (in  future  moist  versions  of  the  model),  its  presence  is 
felt  by  both  scattered  and  absorbable  components.  If  the  cloud  is  thick  enough,  and 
covers  enough  sky,  incident  radiation  can  be  reduced  to  zero.  The  model  allows  for 
variable  amounts  of  cloud  in  each  atmospheric  layer  expressed  as  a  percentage.  Low 
and  middle  level  clouds  can  be  included  as  a  single  layer  covering  all  or  part  of  the 
sky.  The  cloud  layer  directly  affects  both  scattered  and  absorbable  shortwave 
radiation. 


Fractional  absorption  by  water  vapor  is  calculated  by 


ABS(Jc)  =  (0.271) 


4k) 


t0.303 


cos(ZT} 


(3.2) 


where  EH?Q  is  the  effective  amount  of  water  vapor  in  layer  k.  The  radiation  that 
is  finally  absorbed  in  the  soil  becomes  one  component  of  the  surface  energy  balance. 
The  absorbed  part  of  the  incident  radiation  at  the  ground  is 


GLWo  =  (0.349)S0cos(ZT)  -  £ABS(fc)  (3.3) 

k 

where  the  sum  is  taken  over  all  layers.  The  total  absorption  at  the  ground  is  then 


GAB  =  (1  -as)(GLW,  +  GLWo' 


(3.4) 
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To  find  the  IR  flux,  the  equation  of  radiative  transfer  is  solved  subject  to 
the  boundary  conditions  that  downward  IR  flux  at  the  top  of  the  atmosphere  is 
zero,  and  the  upward  IR  flux  at  the  earth’s  surface  is  the  black-body  radiation  at 
the  surface  temperature.  Weighted  transmission  functions  are  used,  corrected  for 
the  pressure  dependence  of  absorption.  The  total  transmission  function  is  assumed 
to  be  the  product  of  the  individual  ones  for  C02  and  H20.  Downward  flux  is 


lRd  *=■  irfii  —  irBoT(iiI  —u*,  Ta)  —  (irB-  —  itBo)T(u*  —  u* ,  T) 


+ 


/xB- 
*BZ 


t (u*-uT,  T)  d(irB) 


(3.5) 


4 

where  for  each  layer,  xB,  =  oT,  .  The  weak  region  is  210-320  K  for  water  vapor, 
so  letting  Tc  =  220  K,  the  weak  dependence  region  need  only  have  a  mean 
temperature  specified  (T).  Similarly,  the  upward  flux  is 


IRu 


XBz  -f- 


T(U*-U ?,  T)  d(xB) 


(3.6) 


and  the  net  upward  flux  is 


IRZ  =  IRu  -  IRd 


(3.7) 


The  only  difficulty  is  determining  the  proper  transmission  function  near  the 
particular  level  where  r  varies  exponentially.  The  model  uses  an  interpolation 
factor  that  is  an  empirical  function  of  pressure,  mixing  ratio  and  layer  thickness. 
This  allows  proper  calculation  of  r  without  a  fine  vertical  mesh.  The  mean 
transmission  functions  are  defined  by  empirical  formulae  at  Tc  =  220  k  and  T  — 
260  K.  Temperature  dependence  of  r  for  CO:  is  neglected,  so  a  mean  r  for  CCk  is 
used  based  on  pressure  and  amount  of  C02.  The  distribution  of  CO.  at  each  level 
is  a  constant. 

The  IR  flux  is  computed  only  for  the  surface,  since  the  IR  cooling  rates  in 
the  free  atmosphere  are  insignificant  on  diurnal  time  scales.  This  saves 
considerable  computation  time.  The  net  radiation  is  then 

NR  =  GAB  -  IR;  (3.8) 
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Radiation  calculations  were  made  with  the  above  scheme  for  a  varying 
number  of  atmospheric  layers.  These  calculations  showed  that  when  layers  were 
thicker  than  100  mb,  more  significant  errors  occurred  in  the  net  radiation  values. 
Comparisons  were  made,  for  example,  using  a  sounding  that  originally  had  19 
unevenly  spaced  levels.  Reducing  this  to  only  7  levels  produced  NR  values  in  error 
by  nearly  15%  compared  to  the  19-level  result.  This  error  is  probably  acceptable 
for  most  situations  since  it  implies  a  fairly  small  error  in  the  surface  temperature 
change  due  to  radiation  and  because  the  model  will  not  be  integrated  for  more  than 
about  24  h.  This  error  can  be  improved  somewhat,  however,  without  the  need  for 
more  detail  in  the  sounding.  If  the  routine  started  with  the  same  7  levels,  but  was 
allowed  to  linearly  interpolate  a  new  level  in  the  middle  of  any  layer  thicker  than 
100  mb,  the  error  in  NR  was  reduced  to  less  than  10% — despite  the  fact  that  no 
additional  vertical  resolution  in  temperature  was  available. 

b.  Surface  energy  balance 

The  surface  energy  balance  has  the  form 

NR  =  SH  +  LH  +  GS  (3.9) 

where  SH  is  the  sensible  heat  flux  upward  from  the  surface,  LH  is  the  latent  heat 
flux  upward  from  the  surface,  and  GS  is  the  soil  heat  flux  downward  into  the 
ground  (heating  the  soil).  The  quantity  NR  is  determined  from  the  radiation 
parameterization  [see  Eq.  (3.8)1. 

The  sensible  heat  flux  and  latent  heat  flux  (SH,  LH)  are  parameterized  using 
Monin-Obukhov  similarity  theory  for  the  planetary  boundary  layer  (PBL).  The 
fluxes  depend  on  the  vertical  gradients  of  temperature  and  specific  humidity  in  the 
surface  layer,  the  wind  in  the  model  layer  directly  above  the  PBL,  and  the  stability 
of  the  boundary  layer.  The  theory  assumes  that  the  structure  of  temperature  and 
moisture  in  the  PBL  have  forms  that  can  be  described  by  universal  structure 
functions  when  scaled  equations  are  used.  There  are  actually  two  structures 
involved,  since  the  PBL  contains  at  least  two  distinct  layers — the  surface  layer  and 
the  boundary  layer  (see  Fig.  3).  Although  the  parameterization  does  not  include  an 
explicit  surface  layer,  one  is  assumed  to  be  present.  This  “surface  layer"  is 
assumed  to  be  a  constant  5  mb  thick.  If  the  functions  are  required  to  be  matched 
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at  their  common  boundary,  the  following  results  are  obtained 


UNSTABLE 


STABLE 


t 

5  mb 

4 


Fig.  3.  Schematic  of  boundary  layer  structure  for  (a)  stable 
and  (b)  unstable  cases. 

The  inverted  equations  have  the  form 

=  CARo,  S )  (3.19) 

PV a‘ 


(3.20) 


C\ tRo,  S, 


(3.21) 


where  Cy(Ko,  S),  C»»(Rc,  S),  and  a0!fio,  S)  have  different  forms  for  stable  and 
unstable  PBLs.  The  latent  heat  flux  (LH)  is  computed  by  assuming  that  it  obeys 
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the  universal  function  for  SH,  that  is, 


Q”  — Qo  _  6**  — 8q 

q*  ~  e * 


where 


<7* 


-LH 

u* 


(3.22) 


(3.23) 


c.  Frictional  drag 

The  frictional  effects  of  the  boundary  layer  are  simulated  with  a  simple 
aerodynamic  drag  formulation  based  on  Sutton  (1953).  In  the  lowest  model  layer,  the 
frictional  term  in  (2.2)  is  given  by 


g  9T  ;ry 
rH  377 


where  the  surface  stress  vector  rav  is  given  by 


(3.24) 


Tiy  =  pC d[V qUI  V0vj)  (3.25) 

with  Cd  ==  u+2 /V 0‘,  and  where  V0  represents  the  magnitude  of  the  wind  in  the 
model  layer  just  above  the  boundary  layer. 


d.  Ground  variables 

Ground  temperature  (TG)  and  ground  wetness  (GW)  are  parameterized  by 
“force  restore”  methods  from  Bhumralker  (1975)  and  Deardorff  (1977).  Following 
Bhumralker  (1975),  heat  conduction  in  the  soil  is  described  by 

3T  glz,t)  KdTg(Z,t) 

3?  *  C  g-2  (3-6) 

where  Tjtz.O  is  the  soil  temperature  at  depth  z.  and  time  t.  We  assume  that  TG  is 
described  by 


TG  *=  T  -f  AT0siniu;() 


(3.27] 
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where  T  is  the  average  temperature  of  the  soil,  and  A Tc  is  the  amplitude  of  the 
variance.  Then  the  solution  of  (3.26)  is 


Tg(z,t )  =  T  4-  AT0exp(  — z /d)[sin(tJt  — z /d)) 


(3.28) 


where  d  =  {2K/cu))l/‘  is  the  depth  at  which  the  amplitude  of  A T0  is  negligible.  For 
an  infinitely  thin  soil  layer,  the  heat  flux  into  the  soil  at  depth  z  is 


G(z,t)  = 


.3  Tg(z,t) 


(3.29) 


Combining  (3.28)  and  (3.29)  gives 


G(z,t )  =  (^^]^AT0  e  Z  ^[sin(u)t  —  z  /  d)  +  cos  (ujt—z/d)} 


(3.30) 


Eliminating  A T0  one  obtains 


G(z,t)  =  a7^|(~)  +  Tg(z,t)  - 


(3.31) 


Consider  a  layer  of  soil  from  the  surface  (z  =  0)  to  some  depth  z.  The  time 
rate  of  temperature  change  for  this  layer  is  given  by 


dTg(Z,t) 


rC(z,f)-GS| 


(3.32) 


If  the  approximation  is  made  that 


Tg(z  —  1  cm,  t)  =»  TG 


(3.33) 


then  (3.32)  becomes  (with  the  use  of  (3.9)), 

c  4-  =  GS  -  (TG  -  T) 


(3.34) 


GS  -  [~f  (TG  -  T)Jc  4-  (g)5 
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(3.35) 


The  surface  soil  moisture  (GW)  is  found  by  assuming  that  it  changes  due  to 
three  main  processes — precipitation,  evaporation,  and  flux  from  below.  The  bulk 
soil  moisture  (GWB)  is  assumed  to  be  constant  over  the  period.  According  to 
Deardorff  (1977)  the  bulk  soil  moisture  changes  over  a  time  scale  of  a  few  weeks,  so 
GWB  can  be  assumed  constant  for  a  24  h  period  with  little  loss  of  accuracy.  The 
surface  soil  moisture  is  changed  according  to 

3GW  _  _  Ci(LH/\-Pr)  _  c2(GW-GWB) 

at  Tr  •  ' 

Deardorff ’s  values  for  the  nondimensional  constants  c.  and  c2  were  computed  from 
data  of  Jackson  (1973),  which  used  measurements  taken  over  bare  soil  near  Phoenix, 
Arizona  in  March.  This  gives 


0.5  GW  ^  75% 

c,  =  14  —  22.5(GW  —0.15)  15%  <  GW  <  75% 

14  GW  £  15% 


c2  =  0.9 


(3.37) 


The  above  discussion  of  the  surface  temperaure  and  soil  moisture 
parameterizations  applies  for  grid  boxes  composed  of  land.  For  grid  points  over 
water,  the  value  of  TG  is  set  to  be  the  temperature  of  the  water  and  held  constant 
over  the  integration,  and  the  soil  moisture  is  set  to  reflect  a  water  surface  capable 
of  continual  evaporation.  Grid  boxes  composed  of  a  mixture  of  land  and  water 
surfaces  should  produce  fluxes  that  represent  some  sort  of  weighted  average 
between  those  from  all-land  boxes  and  those  from  ail-water  boxes.  Several  schemes 
using  the  framework  of  the  above  parameterizations  were  tested  to  try  to  make  use 
of  knowledge  of  the  land-water  mix  of  the  grid  boxes,  but  none  of  them  proved 
satisfactory  for  all  possible  types  of  conditions.  While  a  weighted-average  flux 
approach  seems  very  straightforward,  the  coupling  of  the  surface  temperature 
prediction  with  the  incoming  and  outgoing  radiation  calculation  and  the  forecast 
fluxes  leads  to  approximations  that  cannot  be  good  for  both  day  and  night.  Other 
approaches  led  to  similar  difficulties. 
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e.  Boundary-layer  height 


1)  Unstable  boundary  layer 

The  unstable  PBL  is  assumed  to  be  well  mixed  below  an  inversion 
characterized  by  a  jump  discontinuity  in  potential  temperature  (A0).  The  depth  of 
the  unstable  boundary  layer,  h,  and  the  strength  of  the  inversion,  AS,  are  predicted 
according  to  Zeman  and  Tennekes  (1977).  Their  method  assumes  that  the  PBL  depth 
changes  due  to  turbulent  entrainment  of  air  above  the  inversion  into  the  PBL.  The 
energy  comes  from  the  virtual  SH  flux  at  the  surface,  and  the  change  of  depth  with 
time  depends  on  the  strength  of  the  inversion.  They  use  the  turbulent  kinetic 
energy  (TKE)  budget  to  develop  a  simple  set  of  equations  to  describe  this  process. 

The  sensible  heat  flux  at  the  inversion  is  equal  to  the  temperature  jump,  AS, 
times  the  rate  of  rise  of  the  inversion 


VSH*  =  A0|J  (3.38) 

where  VSH*  is  the  virtual  sensible  heat  flux  at  the  inversion.  The  inversion 
strength  changes  as  a  function  of  entrainment  of  stable  air  from  above,  and  net 
sensible  heat  transfer  inside  the  boundary  layer.  It  is  given  by 


=  7  If  “  (VSH  - 


(3.39) 


where  VSH  is  the  virtual  sensible  heat  flux  at  the  surface  and  7  is  the  potential 
temperature  lapse  rate  above  the  PBL.  The  TKE  budget  can  be  written  as 


3(TKE) 
3 1 


—  production  +  transfer  -f  dissipation 


(3.40) 


which  can  be  expanded  into 


C  ^  *  3h  _  Q  v^h  i  c ~  ^ * 

w  —  ^7  —  vili*  ■+  L.f 


h  3 1 


TS 


—  C d  u.*  u>t 


(3.41) 


where  Ct,  C and  Cd  are  dimensionless  constants.  Substituting  for  dh /dt  from 
(3.38)  yields 
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-  TS  VSH-  "  X  [C/-Cau,tJX)] 

Substituting  for  w+  (see  appendix  A)  gives 

vsHh  _  r  r  u .  \  h 

— vsFT  _ 

The  values  of  the  dimensionless  coefficients  Cd,  Cy,  and  Ct  are  taken  to  be  (Zeman 
1975) 


14 


Ctu>iTS 


gA6h 


(3.43) 


Ctu>»  TS 
gA&h 


(3.42) 


Cd  •=  0.024 

Cy  =  0.50  (3.44) 

Ct  =  3.55 


We  can  write  the  rate  of  change  of  boundary  layer  height  as 

dh  =  _  VSH* 

3 1  A0 


(3.45) 


In  the  case  where  A0  =  0,  no  inversion  exists  and  the  atmosphere  presents  no 
barrier  to  inversion  rise.  In  this  case,  the  model  assumes  a  very  small  value  for 
A0,  since  the  inversion  must  rise  at  a  rapid  but  finite  rate  due  to  turbulent 
entrainment. 


2)  Stable  boundary  layer 

The  depth  of  the  stable  boundary  layer  is  calculated  using  a 
parameterization  from  Zeman  (1979).  The  mean  kinetic  energy  budget  is  the  basis 
for  a  rate  equation  governing  the  height  of  the  stable  boundary  layer.  This  results 
in 


+  ^  (3.46) 

|Au  '  dr 

where  the  quantity  CP*  is  a  tunable  function  of  the  Richardson  number,  determined 
by  comparison  with  numerical  and  laboratory  results,  and  (  is  the  angle  between  the 
geostrophic  wind  above  the  boundary  layer  and  the  surface  stress. 


dh 

dt 


C; r 


il  4-5 h  L) 


3|Aur 


(Aucosf  4-  At’sinf) 
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/.  Boundary-layer  variables 


The  surface  temperature  is  found  as  a  by-product  of  the  PBL  depth 
calculations.  For  the  unstable  PBL,  we  proceed  as  follows: 

A 0P  =  e,  +  Zt~A/l  (flp-fl,)  -  flP  (3.47) 

where  A h  is  the  change  in  PBL  depth  over  one  time  step. 

The  potential  temperature  at  the  top  of  the  5  mb  surface  layer  is  given  by 


es  =  ep  -  a 9 


(3.48) 


which  allows  the  surface  temperature  to  be  found  by  Poisson’s  equation  as 


TS  =  0= 


Ps/c  5 


tO.286 


1000 


(3.49) 


For  the  stable  PBL,  we  start  in  a  similar  manner  pt  that  there  is  no 

inversion  jump  discontinuity.  We  simply  interpolate  in  potential  temperature 
between  the  PBL  top  and  the  ground  to  find  then  compute  TS  with  (3.49). 


g.  Boundary -layer  transitions 

Transitions  between  stable  and  unstable  boundary-layer  regimes  have  been 
incorporated  into  the  boundary-layer  parameterization.  In  the  real  atmosphere, 
these  transitions  are  almost  certainly  abrupt  and  characterized  by  poorly  defined 
structures.  Within  the  model,  such  a  situation  would  lead  to  numerical  problems. 
In  addition,  there  is  no  known  parameterization  for  a  boundary  layer  in  transition, 
so  we  have  to  provide  a  mechanism  for  an  orderly  transition. 

The  transition  stage  is  determined  by  the  sign  of  the  sensible  heat  flux. 
When  the  direction  of  the  sensible  heat  flux  is  incompatible  with  the  nature  of  the 
boundary  layer,  the  transition  stage  is  set.  The  boundary-layer  height  is 
constrained  to  fall  rapidly  from  its  current  location  to  a  height  of  90  m  (arbitrarily 
chosen).  To  do  this,  the  boundary-layer  height  tendency  is  initially  set  at  a  value 
that  would  cause  the  boundary  layer  to  be  90  m  after  15  min.  The  tendency  is 
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adjusted  each  time  step  according  to 


dh 

dt 


dhl  fh  -  90p 
[dt\old{  90 


—  (6  m  min  ') 


(3.50) 


where  the  subscript  “old”  refers  to  the  previous  timestep.  The  multiplicative 
factor  slows  the  rate  of  change  as  the  transition  nears  completion  while  the  6  m 
min-1  additive  constant  assures  that  the  transition  will  be  completed  in  a  finite 
time.  During  this  transition,  the  potential  temperature  structure  in  the  boundary 
layer  is  changed  as  little  as  possible,  thereby  allowing  the  proper  parameterization 
to  begin  with  the  current  structure  as  though  the  transition  had  been 
instantaneous.  The  potential  temperature  at  the  top  of  the  boundary  layer  is 
linearly  interpolated  between  its  current  value  and  the  potential  temperature  of  the 
ground.  Boundary-layer  moisture  is  conserved  during  this  process  of  boundary- 
layer  height  fall,  again  to  preserve  current  conditions  as  much  as  possible.  When 
the  boundary-layer  height  reaches  90  m  the  transition  ends.  The  boundary  layer 
then  begins  to  grow  again  according  to  the  stability  at  that  time. 

Transitions  from  unstable  to  stable  and  vice  versa  have  been  tested  and 
found  to  be  relatively  smooth  and  reliable.  Although  most  of  these  transitions  will 
take  place  at  dawn  or  dusk,  it  is  conceivable  that  a  major  change  in  cloud  cover  or 
air  mass  might  initiate  such  a  change  at  any  time  during  an  integration.  The  only 
aspect  of  this  that  will  not  be  handled  well  by  the  parameterization  is  a  situation  in 
which  a  well-mixed  convective  boundary  layer  grows  to  a  substantial  depth,  then 
the  radiation  is  interrupted  by  clouds,  causing  a  transition  to  a  stable  boundary 
layer,  followed  by  a  breakup  of  the  clouds  and  a  return  to  unstable  stratification. 
The  new  convective  boundary  layer  will  start  over  and  the  atmosphere  above  the 
boundary  layer,  which  in  the  real  atmosphere  would  be  almost  dry  adiabatic  in 
structure,  will  be  slightly  stable  in  the  model.  This  is,  of  course,  a  problem  with 
the  limited  vertical  resolution  of  the  model,  not  the  fault  of  the  parameterization. 
Only  extensive  testing  will  be  able  to  show  if  this  is  a  serious  problem.  We  suspect 
that  most  cloud  cover  will  only  reduce  the  sensible  heat  flux  in  magnitude  rather 
than  actually  reversing  the  sign.  Hence,  the  boundary  layer  will  grow  much  more 
slowly  but  will  not  go  through  a  transition. 
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4.  Simulations  with  the  model 


a.  Mountain  lee-wave  simulations 

A  critical  test  of  any  meso-0  scale  model  is  the  simulation  of  mountain  lee- 
waves  (Anthes  and  Warner  1978;  Nickerson  et  al.  1986).  While  previous 
investigators  have  usually  made  these  tests  with  a  two-dimensional  version  of  their 
models  (Anthes  and  Warner  1978;  Nickerson  et  al.  1986),  we  choose  to  use  the  full 
three-dimensional  model — except  without  any  boundary  layer  or  radiation 
parameterizations.  By  choosing  a  long  uniform  ridge  oriented  perpendicular  to  the 
flow,  we  can  obtain  a  nearly  two-dimensional  solution  crossing  the  center  of  the 
ridge  (which  allows  comparison  with  previous  two-dimensional  simulations),  while  not 
modifying  the  three-dimensional  staggering  of  the  variables  as  required  in  a  two- 
dimensional  analog  model.  Since  the  boundary-layer  parameterization  is  turned  off 
for  this  test,  the  boundary-layer  height  must  be  specified  and  is  held  constant. 
While  many  boundary-layer  heights  were  tested,  including  artificially  varying  ones, 
we  will  show  the  results  when  crh  =  0.80  which  represents  equally  spaced  layers  for 
the  model. 

For  the  simulation  discussed  in  this  section,  the  thermodynamic  profile  was 
specified  as  the  U.  S.  Standard  Atmosphere,  the  specific  humidity  was  set  to  zero, 
and  the  u-component  winds  were  specified  as  20  m  s-1  at  all  levels.  The  Coriolis 
force  is  included  in  the  model  for  these  simulations,  and  the  initial  conditions  were 
in  approximate  geostrophic  balance  (at  a  latitude  of  about  40"N).  A  mountain  ridge 
with  Gaussian  cross  section  and  a  maximum  elevation  of  1  km  was  oriented 
north-south  near  the  center  of  the  FGM  domain.  The  ridge  had  a  constant  height 
for  the  center  12  north-south  grid  points  and  decreased  linearly  to  zero  height  so 
that  it  did  not  extend  into  the  CGM  collar.  This  avoided  problems  associated  with 
the  specification  of  the  solution  on  the  CGM  lateral  boundary.  Though  not 
presented  below,  simulations  with  the  ridge  rotated  90"  and  with  winds  from  the 
south,  north,  and  west  were  carried  out  to  isolate  any  coding  errors. 

Figure  4  shows  a  vertical  cross  section  through  the  FGM  domain  after  12  h 
simulated  time.  The  lsentropic  cross  section  shows  the  expected  lee-wave  structure 
with  the  waves  sloping  upstream  with  height;  this  is  also  reflected  in  the  horizontal 
velocity  field.  While  the  isentropes  in  Fig.  4  are  plotted  objectively  by  linear 
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interpolation,  the  isotachs  are  subjectively  analyzed  because  of  the  coarse  vertical 
resolution.  These  results  are  quite  similar  to  those  obtained  for  similar  conditions 
by  Nickerson  et  al.  (1986)  using  a  model  with  16  layers,  and  by  Anthes  and  Warner 
(1978)  using  a  model  with  only  6  layers.  It  should  be  noted  that  Anthes  and  Warner 
(1978)  discussed  the  requirement  of  a  relatively  large  horizontal  diffusion  in  order 
for  a  model  with  few  layers  and  no  explicit  damping  region  near  the  upper  boundary 
to  represent  lee-waves  successfully.  As  discussed  above,  the  value  of  K  used  here 
leads  to  a  nondimensional  diffusion  parameter  that  meets  the  criteria  established  by 
Anthes  and  Warner  (1978). 


Fig-  4.  Vertical  cross  section  throught  the  FGM  domain  at  12  h  simulated  time  for 
lee  wave  simulation  showing  the  potential  temperature  (solid),  and  the  u-component 
velocity  (dashed),  Isentropes  are  labeled  in  K  along  the  right  and  the  velocity  is 
contoured  in  m  s-1. 
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Fig.  5.  Horizontal  plot  for  simulation  shown  in  Fig.  3.  Solid  contours  give 
sea-level  pressure  (in  mb),  thin  dashed  contours  are  for  surface  topography 
(in  m)  of  the  mountain  ridge,  and  boundary  layer  winds  (one  full  barb 
represents  5  m  s~')  are  plotted  at  each  grid  point. 
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Figure  5  shows  a  plot  of  the  surface  pressure  field  and  boundary-layer  winds 
for  this  simulation  at  12  h  simulated  time.  This  figure  shows  a  lee  trough  and 
upstream  flow  blocking  that  is  often  associated  with  flow  over  a  ridge.  It  is  quite 
evident  that  the  model  solution  is  smooth  and  well-behaved,  even  at  the  lateral 
boundaries  representing  the  FGM/CGM  interface.  Figure  6  shows  the  average  rate 
of  change  of  surface  pressure  in  the  FGM  during  the  integration.  It  is  clear  that 
the  model  undergoes  an  adjustment  during  the  first  few  hours  of  integration  as  the 
pressure  adjusts  to  the  mountain  lee-wave  structure.  The  model  is  quite  steady  for 
the  last  6  hours  of  the  integration,  a  fact  confirmed  by  plots  (not  shown)  of  winds 
and  surface  pressure  at  8  h  and  10  h. 


Fig.  6.  Average  rate  of  change  of  ir  versus  time  for  FGM  domain  during 
the  simulation  shown  in  Figs.  4  and  5. 
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Although  written  in  the  framework  of  boundary-layer  coordinates,  the  model 
reduces  to  a  traditional  cr-coordinate  model  if  the  boundary-layer  height,  oh,  is 
specified  as  a  constant  in  the  domain — as  it  was  in  the  simulation  discussed  above. 
Numerous  simulations  in  which  the  boundary-layer  height  was  artificially  raised  or 
lowered  were  carried  out  for  the  mountain  lee- wave  case.  Even  when  the  boundary- 
layer  height  was  changing  relatively  rapidly,  the  lee-wave  structure  was  preserved 
and  the  model  maintained  a  near  geostrophic  balance  for  pressure  and  wind  fields 
(Colby  and  Seitter  1990). 

b.  Sea-breeze  simulation 

The  full  capabilities  of  the  boundary-layer  coordinate  system  can  only  be 
tested  if  a  mesoscale  circulation  that  depends  on  boundary-layer  processes  is 
simulated.  A  natural  candidate  for  such  a  test  is  a  sea-breeze  circulation.  In  order 
to  isolate  boundary-layer  processes  from  topographic  forcing,  we  choose  to  simulate 
the  sea-breeze  circulation  of  southern  Florida,  which  we  treat  as  a  flat  land  mass  at 
zero  elevation.  We  will  show  here  a  simulation  similar  to  the  “southeast  flow”  case 
treated  by  Pielke  ( 1 972)-  which  matched  his  initial  conditions  as  much  as  possible. 
This  included  specifying  a  relatively  stable  environment  [using  the  Peilke  (1972) 
sounding]  and  a  7  m  s~‘  synoptic  wind  from  135'  for  all  levels  that  was  in 
approximate  geostrophic  balance  with  the  initial  pressure  field.  The  FGM  domain 
with  the  initialized  boundary-layer  winds  is  shown  in  Fig.  7.  The  head  of  each  wind 
vector  represents  the  center  of  a  thermodynamic  gridbox  in  the  FGM.  There  are 
four  gridpoints  in  Lake  Okeechobee  which  are  specifed  as  water  points  in  the  model. 
The  model  simulation  is  started  at  0500  local  time  on  Julian  day  180  (01  June),  so 
the  simulation  begins  abo  one  hour  before  sunrise.  Soil  moisture  content  is  set 
uniformly  on  all  land  grid  points  as  50%,  and  the  surface  roughness  length  is  set  on 
land  and  water  as  0.04  m  and  0.005  m,  respectively.  The  initial  boundary-layer 
height  is  specified  as  ah  —  0.96  over  the  entire  FGM  and  CGM  domain,  which 
translates  to  a  boundary-layer  depth  of  about  320  m.  This  simulation  is  not  meant 
to  reproduce  a  specific  observation  of  the  Florida  sea-breeze,  but  the  conditions  are 
representative  of  a  class  of  synoptic  situations  that  leads  to  an  inland  convergence 
zone  over  the  Florida  penninsula  (Peilke  1972). 

Figures  8-10  show  the  boundary-layer  winds  in  the  FGM  domain  at  6  h,  8  h, 
and  10  h.  The  major  features  of  the  mesoscale  structure  are  quite  similar  to  those 
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Fig.  7.  FGM  domain  and  initial  boundary-layer  winds  for  Florida  sea- 
breeze  simulation. 
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Fig.  9.  FGM  boundary  layer  winds  at  8  h,  for  the  Florida  sea-breeze 
simulation . 
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Fig.  10.  FGM  boundary  layer  winds  at  10  h,  for  the  Florida  sea-breeze 
simulation. 
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presented  by  Pielke  (1972).  The  sea-breeze  initiated  on  the  western  coast  must 
work  against  the  synoptic  scale  flow  while  the  sea-breeze  circulation  on  the  eastern 
coast  simply  augments  the  southeasterly  flow  and  turns  it  somewhat.  A 
convergence  line  is  established  near  the  western  coast.  While  the  20  km  resolution 
of  the  present  model  is  not  capable  of  reproducing  the  complicated  wake-like 
structure  downstream  of  Lake  Okeechobee  simulated  by  Pielke’s  10  km  resolution 
model,  there  is  clear  modification  of  the  flow  by  the  weaker  lake-breeze  in  that 
region. 

An  east-west  vertical  cross-section  through  the  model  at  9  h  simulated  time 
is  shown  in  Fig.  11.  This  cross-section  was  taken  at  grid  row  15,  which  passes 
through  the  center  of  Lake  Okeechobee.  The  boundary-layer  height,  denoted  by 
the  thin  dashed  line,  is  quite  high  over  the  warm  land  regions  due  to  growth  as  an 
unstable  boundary  layer.  The  boundary  layer  remains  stable  and  low  over  the 
water  surfaces,  which  are  much  cooler  than  the  land  at  this  time.  The  temperature 
gradient  near  the  eastern  coast  is  weaker  than  that  near  the  western  coast  because 
cooler  air  is  advecting  in  from  the  east  and  being  warmed  by  surface  heat  fluxes 
from  the  land.  At  the  western  coast,  where  the  sea-breeze  circulation  is  just  able 
to  overcome  the  synoptic  winds,  there  is  little  advection  so  a  large  gradient  is 
established  between  the  cool  marine  air  and  the  warm  air  over  the  land.  The 
structure  of  the  boundary  layer  in  the  vicinity  of  the  coast  is  very  similar  to  that 
produced  by  Anthes  (1978)  with  a  model  utilizing  many  layers  near  the  surface  in 
order  to  resolve  explicitly  the  boundary  layer  changes  during  the  development  of 
the  sea-breeze.  Despite  the  large  gradients  in  the  boundary-layer  height  near  the 
coast,  very  little  noise  is  generated.  This  is  a  result  cf  the  well-posedness  of  the 
boundary-layer  coordinate  system  and  the  horizontal  diffusion  operating  on  the 
variables. 

Figure  12  shows  the  boundary  layer  winds  in  the  FGM  at  24  h.  After  sunset, 
the  boundary  layer  undergoes  transition  to  a  stable  boundary  layer,  and  as  the 
temperature  difference  between  the  land  and  water  decreases,  the  winds  over  the 
whole  domain  tend  to  evolve  back  to  the  initial  southeasterly  flow.  The  synoptic 
flow  over  the  ocean  has,  by  this  time,  turned  to  be  somewhat  more  easterly.  This 
turning  toward  lower  pressure  occurs  over  the  course  of  the  integration  as  the 
boundary-layer  winds  over  the  ocean  establish  a  new  balance  from  the  geo6trophic 
initial  conditions  to  one  that  includes  friction.  This  is  not  evident  in  the  northern 
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Fig.  12.  FGM  boundary  layer  winds  at  24  h,  for  the  Florida  sea-breeze 
simulation. 
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portion  of  the  FGM  domain  over  the  oceans  because  a  substantial  portion  of  the 
northern  CGM  collar  included  land  and  the  sea-breeze  circulation  established  there 
had  a  large  impact  or.  the  CGM  influence  on  the  FGM  winds.  The  lend  has  cooled 
radiatively  during  the  night  hours  and  at  this  time  is  slightly  cooler  than  the  ocean 
surface — leading  to  a  slight  land  breeze  evident  in  Fig.  12  as  reduced  velocity  flow 
on  the  eastern  coastal  region.  There  is  also  some  evidence  in  the  model  results  to 
suggest  that  the  reduction  of  easterly  component  in  the  eastern  part  of  the  domain 
is  a  result  of  eastward  momentum  that  remained  in  the  boundary-layer  flow  from 
the  west  coast  sea-breeze  and  propagated  across  the  penmnsula.  The  stable 
boundary  layer  over  land  at  this  time  has  a  height  of  about  160  m. 


S.  Initial  tests  using  realistic  complex  terrain 

a.  Test  domain 

Future  research  will  involve  the  testing  of  a  prototype  model  using  real  data 
in  clear  weather  and  stratiform  precipitation  situations.  Such  testing  requires  a 
specific  geographic  location.  Ideally,  this  location  should  experience  a  wide  variety 
f  weather  phenomena  such  as  sea-breeze  circulations,  synoptic  scale  stratiform 
precipitations  events,  orographic  enhancement  of  precipitation,  and  other  mesoscale 
patterns  for  which  the  model  is  intended  to  provide  guidance.  Southern  New 
England  represents  one  of  several  locations  that  could  be  considered,  and  given  our 
obvious  desire  for  operational  guidance  in  our  own  area,  we  have  chosen  this 
geographic  region  for  testing  of  the  prototype  model. 

The  domain  for  the  fine  grid  mesh  is  shown  in  Fig.  13  (only  the  locations  of 
“thermodynamic”  grid  points  are  shown).  The  coarse  grid  mesh  which  surrounds  the 
FGM  domain  is  shown  in  Fig.  14.  In  addition  to  being  a  coastal  location,  this  area 
has  complex  terrain  that  is  rich  in  detail  on  meso-0  scales.  This  is  quite  evident  in 
the  contour  plot  of  the  FGM  topography  shown  in  Fig.  15.  The  procedure  followed 
to  create  the  partial-envelope  orography  displayed  in  this  figure  is  described  in 
detail  in  appendix  C.  Several  mesoscale  orographic  features  thought  to  be 
important  to  the  weather  in  the  New  England  area  are  obvious.  These  include  the 
Adirondack  mountains,  the  Hudson  River  Valley,  the  Green  Mountains,  the 
Berkshires,  the  Connecticut  River  Valley,  the  White  Mountains,  and  the  Worcester 
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FINE  9RIQ  MESH  OOHRIN 

Fig.  13.  FGM  domain  for  New  England  test  simulations  showing  the 
location  of  thermodynamic  points. 
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Fig.  14.  CGM  domain  for  New  England  test  simulations 
location  of  thermodynamic  points. 


SURTflCE  CLcVflTION  IH1 


Fig.  15.  FGM  terrain  field  (see  appendix  C  for  details  on  the  creation  of 
this  partial  envelope  terrain  field). 
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Fig.  16.  Map  of  FGM  domain  showing  major  geographic  features. 


41 


Plateau  (see  Fig.  16  for  a  map  showing  the  locations  of  these  features).  A  potential 
problem  is  that  Mt.  Washington  falls  almost  directly  on  a  grid  point  at  the  northern 
edge  of  the  FGM  domain  and  therefore  tends  to  dominate  the  terrain  field  in 
northern  New  Hampshire.  This  particular  point  may  need  to  be  reduced  in  height  in 
order  to  avoid  numerical  problems  near  the  FGM  boundary. 

As  discussed  in  section  3,  the  model’s  boundary-layer  parameterization  scheme 
requires  a  grid  box  to  be  composed  entirely  of  either  land  or  water,  but  the  surface 
roughness  parameter  can  be  adjusted  to  be  consistent  with  an  average  value  for  a 
grid  box  with  a  mixture  of  land  and  water.  This  means  some  of  the  surface 
roughness  effects  of  the  smaller  scale  features  can  be  included  in  the  model,  even 
though  the  coastline  is  a  “blocky”  approximation  to  the  actual  coastline  shown  in 
Fig.  13  for  use  in  the  surface  flux  calculations.  U.S.  Geologic  Survey  maps  were 
used  to  make  a  careful  subjective  estimate  of  the  percent  of  water  coverage  for 
each  grid  box,  centered  on  the  thermodynamic  grid  points  shown  in  Fig.  13.  A 
contour  plot  of  percent  water  coverage  for  the  FGM  domain  that  resulted  from  this 
procedure  is  shown  in  Fig.  17.  Notice  how  well  even  subtle  variations  in  the 
coastline  are  captured  in  the  contouring,  and  the  influence  of  some  important  inland 
bodies  of  water,  such  as  the  Hudson  River,  the  Quabbin  Reservoir,  and  Lake 
Winnipesaukee  (see  Fig.  16).  For  the  boundary-layer  parameterization,  a  grid  point 
is  considered  to  be  water  if  its  grid  box  is  composed  of  more  than  50%  water 
coverage.  The  island  of  Nantucket  presented  a  problem  in  this  analysis.  Review  of 
Fig.  13  shows  that  Nantucket  lies  nearly  in  the  middle  of  a  square  formed  by  four 
grid  points.  This  means  that  its  land  is  almost  equally  divided  between  the  four 
grid  boxes  represented  by  those  four  grid  points.  Rather  than  have  the  model 
perceive  a  40X40  km  area  of  partial  land  coverage,  the  decision  was  made  to  “move” 
Nantucket  approximately  10  km  northwest  so  that  its  land  mass  fell  entirely  within 
one  grid  box.  This  move  is  obvious  in  Fig.  17,  in  which  the  partial  land  contour 
representing  Nantucket  is  displaced  northwestward  of  the  geographic  location  of  the 
island. 

The  initial  smoothed,  partial  envelope  terrain  field  for  the  CGM  is  shown  in 
Fig.  18.  During  the  first  test  simulations  using  the  New  England  domain,  it  became 
apparent  that  there  was  an  incompatibility  between  the  CGM  terrain  field  and  the 
flow  relaxation  boundary  condition  used  on  the  CGM  lateral  boundaries.  The 
results  of  these  simulations  and  the  changes  in  the  CGM  terrain  that  resulted  are 
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Percent  water  coverage  in  each  grid  box  for  FGM  domain. 
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Fig.  18.  Unmodified,  smoothed,  partial  envelope  terrain  field  for  the  CGM 
domain. 
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Fig.  19.  Percent  water  coverage  in  each  grid  box  for  CGM  domain 
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discussed  in  the  next  section.  A  fractional  water  coverage  field  for  the  CGM  was 
constructed  in  a  manner  similar  to  that  used  for  the  FGM,  and  this  field  is  shown  in 
Fig.  19. 

b.  Initial  simulations  and  modifications  to  CGM  terrain  field 

A  very  simple  set  of  initial  conditions  was  chosen  for  these  first  test 
simulations  of  the  model  on  realistic  terrain.  The  thermodynamic  and  moisture 
profiles  were  initialized  as  the  U.S.  Standard  Atmosphere  temperature  structure  at 
all  locations  with  a  relative  humidity  of  50%.  The  winds  were  initialized  as 
westerly  at  all  heights  with  a  magnitude  4  m  s-',  and  the  surface  pressure  field  was 
set  to  be  in  near-geostrophic  balance  with  these  winds.  The  simulation  was  started 
at  0530  local  time  on  01  June,  which  is  shortly  after  sunrise  for  the  center  of  the 
domain.  The  boundary-layer  height  was  initialized  at  oh  =  0.96  over  the  entire 
domain,  yeilding  a  boundary-layer  thickness  of  just  over  300  m. 

Figure  20  shows  the  sea  level  pressure  and  boundary  layer  wind  fields  in  the 
CGM  domain  after  1  h  simulated  time.  Comparison  with  Fig.  18  makes  it  clear  that 
large  pressure  anomalies  are  being  created  at  grid  points  near  the  boundary  where 
the  terrain  height  differs  substantially  from  the  terrain  height  of  a  boundary  point. 
The  explanation  for  this  is  really  quite  straightforward.  The  flow  relaxation 
boundary  condition  employs  a  Newtonian  damping-type  term  to  relax  the  interior 
value  of  a  quantity  toward  the  value  of  that  quantity  on  the  nearest  boundary 
point  (Davies  1976).  A  weighting  function  controls  the  amount  of  damping  in  the 
relaxation  region,  smoothly  decreasing  the  influence  of  the  boundary  point  as  the 
distance  from  the  boundary  increases  over  a  5Ax-wide  region.  Therefore,  surface 
pressure  values  in  the  relaxation  region  will  be  forced  toward  the  boundary  value 
of  surface  pressure.  If  the  boundary  point  is  at  a  different  elevation  than  the 
interior  point,  the  forcing  will  tend  to  relax  the  interior  point  to  the  boundary 
surface  pressure  even  though  this  value  is  not  consistent  in  terms  of  sea  level 
reduced  pressure.  While  it  may  appear  that  this  situation  could  be  corrected  using 
sea  level  pressure  in  the  relaxation  condition  rather  than  surface  pressure,  similar 
height  dependencies  exist  for  the  other  variables,  such  as  temperature  and  moisture, 
with  similar  negative  impacts  on  the  simulation. 

Problems  arising  from  use  of  the  relaxation  condition  in  complex  terrain  near 
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Fig.  20.  Sea  level  pressure  and  boundary  layer  wind  fields  for  CGM 
domain  after  1  h  simulation  when  the  terrain  field  shown  in  Fig.  18  is  used. 


the  boundary  can  be  reduced  by  greatly  smoothing  the  terrain  field  in  the 
relaxation  region.  The  only  way  to  truly  eliminate  this  source  of  accumulating 
error  is  to  set  the  terrain  height  of  each  point  in  the  relaxation  region  to  be  equal 
to  the  height  of  its  nearest  boundary  point.  To  do  this,  the  terrain  heights  of  the 
innermost  grid  points  in  the  relaxation  region  (5Ai  from  the  boundary)  were 
extended  outward  to  the  boundary.  This  allows  variation  in  the  terrain  heights 
along  the  boundary  and  produces  the  “corrugated”  outer  terrain  field  shown  in  Fig. 
21.  The  fact  that  this  modified  terrain  field  does  not  match  a  smoothed  actual 
terrain  for  the  outermost  region  of  the  CGM  reinforces  the  notion  that  the  solution 
.n  ♦  relaxation  region  should  not  be  interpreted  as  a  physical  solution,  but  as  a 
modified  solution  that  isolates  the  interior  from  the  negative  impacts  of  the  lateral 
boundaries.  Care  will  need  to  be  taken,  however,  in  specifying  the  large-scale 
forcing  on  the  outer  CGM  lateral  boundaries  since  the  terrain  heights  at  these  grid 
points  will  not  necessarily  match  the  terrain  heights  in  the  dataset  supplying  the 
external  forcing.  This  source  of  error  was  already  present,  of  course,  since  the 
large-scale  forcing  data  will  come  from  a  much  coarser  resolution  grid  whose 
smoothed  terrain  field  would  not  necessarily  match  that  of  the  of  the  CGM  even 
without  the  modification  described  above. 

Figure  22  shows  the  sea  level  pressure  and  boundary  layer  wind  fields  after 
1  h  simulated  time  for  the  same  initial  conditions  used  to  produce  Fig.  20,  but  when 
the  modified  terrain  field  shown  in  Fig.  21  is  used.  The  solution  is  obviously  much 
smoother,  and  the  large  scaie  forcing  that  the  CGM  provides  for  the  FGM  appears 
much  more  natural.  The  FGM  solution  for  this  same  simulation  at  this  same  time  is 
shown  in  Fig.  23.  Even  in  these  early  stages  of  the  simulation,  the  effects  of  the 
terrain  are  visible  in  the  winds  and  pressure  patterns.  The  Mt.  Washington  area 
near  the  northern  boundary  of  the  FGM  appears  to  be  exhibiting  some  perturbations 
in  pressure  that  may  be  a  combination  of  physically  realistic  lee-wave  phenomena 
and  aliasing  at  the  CGM/FGM  interface.  As  mentioned  above,  some  special 
smoothing  may  need  to  be  applied  to  the  terrain  in  this  area  because  of  its 
closenes'-  to  the  edge  of  the  FGM. 
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Fig.  21.  Modified  terrain  field  for  the  CGM  domain  that  sets  the  terrain 
height  of  points  in  the  relaxation  region  the  same  as  the  nearest  boundary 
point  for  consistency  with  the  relaxation  boundary  condition. 


—  49 


G»  VECTOR  WINDS 


CGN  SURERCE  PRESSURE  t*BI 


Fig.  22.  Sea  level  pressure  and  boundary  layer  wind  fields  for  CGM 
domain  after  1  h  simulation  when  the  modified  terrain  field  shown  m  Fig. 
21  is  used. 
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Fig.  23.  Sea  level  pressure  and  boundary  layer  wind  fields  for  FGM 
domain  after  1  h  simulation  when  the  modified  terrain  field  shown  in  Fig. 
21  is  used. 
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6.  Conclusions  and  future  work 


This  paper  has  presented  a  mesoscale  model  developed  by  transforming  the 
equations  of  motion  into  a  coordinate  system  in  which  the  boundary-layer  top  is  a 
coordinate  surface.  This  leads  to  a  means  of  parameterizing  boundary-layer 
processes  in  a  simple  and  computationally  efficient  manner  while  still  incorporating 
those  processes  that  should  be  important  on  the  mesoscale.  Tests  with  the  model 
demonstrate  that  the  boundary-layer  coordinate  system  that  has  been  developed  is 
robust  and  can  model  the  real  atmosphere  reasonably  well,  despite  having  a  very 
limited  vertical  resolution.  In  addition,  the  boundary-layer  parameterizations, 
though  simple,  forecast  the  gross  characteristics  of  both  stable  and  unstable 

boundary  layers  well. 

Additional  work  is  necessary  to  test  this  dry  prototype  model  using  the 
southern  New  England  complex  terrain.  The  ability  of  the  model  to  simulate  the 
interaction  of  orographic  and  sea-breeze  type  phenomena  has  not  been  fully 

evaluated.  Several  nonprecipitating  case  studies  (both  sea-breeze  cases  and  “dry” 
cold-frontal  passages)  have  been  identified  and  will  be  used  to  test  the  model  with 
real  data.  For  these  simulations,  synoptic  scale  forcing  can  be  included  on  the  CGM 
lateral  boundaries  by  interpolating  the  observed  fields  in  both  time  and  space.  In  an 
operational  setting,  synoptic  scale  boundary  forcing  would  be  provided  by  output 
from  one  of  the  operational  synoptic  scale  models. 

While  a  dry  mesoscale  model  such  as  the  one  presented  here  has  many 

operational  uses,  it  is  clear  that  many  of  the  situations  in  which  a  mesoscale  model 
could  provide  added  forecast  guidance  involve  production  and  enhancement  of 

precipitation  by  local  circulations  and  complex  terrain.  Work  is  currently  underway 
to  add  moist  physics  to  the  model  described  here  so  that  it  can  be  used  as  a  more 
general  forecasting  tool. 
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APPENDIX  A 


ABS(fc) 

c 

C 

Cp 

CP* 

Cd,  C f 

Ci,  c2 
C(m) 
C,(Mf) 
di 

Eh2o 

f 

F 

0 

GAB 

GLWc 

GLWS 

GS 

GWB 

GW 

h 

IRd 

IRu 


List  of  Variables 

fractional  absorption  of  radiation  by  water  vapor  for  layer  k 
volumetric  heat  capacity 
dimensionless  constant 

specific  heat  of  dry  air  at  constant  pressure 
a  tunable  function  of  the  Richardson  number 
Ct  dimensionless  constants  taken  from  Zeman  (1975) 
nondimensional  constants 

universal  function  for  a  stable  boundary  layer 

universal  function  for  an  unstable  boundary  layer 

depth  of  diurnal  cycle  (=  10  cm) 

effective  amount  of  water  vapor  in  layer  k 

Coriolis  parameter 

friction  term 

acceleration  of  gravity 

total  absorption  of  radiation  at  the  ground 
absorbed  part  of  incident  solar  radiation 
scattered  part  of  incident  solar  radiation 
soil  heat  flux  downward  into  the  ground 
percent  bulk  soil  saturation  (top  50  cm) 
percent  surface  soil  saturation 
height  of  the  boundary  layer 
downward  flux  if  infra-red  radiation 
upward  flux  if  infra-red  radiation 
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k 


von  Karman’s  constant 


K  eddy  viscosity 

L  Obukhov’s  length 

LH  latent  heat  flux  at  top  of  the  surface  layer 

NR  net  radiation  at  the  surface 

Pr  precipitation  rate 

Q  heating  term 

<7*  -LH /u* 

S0  solar  constant  as  a  function  of  day  of  year 

SH  sensible  heat  flux  at  top  of  the  surface  layer 

SH„  sensible  heat  flux  at  the  boundary-layer  inversion 

T  average  temperature  of  the  soil,  assumed  to  be  invariant  with 

depth  (in  practice,  the  50  cm  soil  temperature  is  used) 

TS  temperature  at  the  top  of  the  surface  layer 

T c  critical  temperature  dividing  the  region  of  weak  temperature 

dependence  of  r  to  that  of  strong  dependence  of  r 

Tg(z,1)  soil  temperature  at  depth  z  and  time  t 

u *  friction  velocity 

u^oi  vg0  components  of  the  surface  geostrophic  wind 

V0  vector  speed  just  above  the  boundary  layer 

VQ  windspeed  just  above  the  boundary  layer 

VSH  virtual  sensible  heat  flux  at  the  surface 

VSH,,  virtual  sensible  heat  flux  at  the  boundary-layer  inversion 

V  vector  velocity  ( =  ui  -+  vj)  , 

u)*  convective  velocity  scale  —  VSH  hj' 

Umax  field  capacity  soil  moisture 

z0  roughness  length  (function  of  location) 

zt  thickness  of  the  next  model  layer  above  the  PBL 

ZT  zenith  angle  for  time  of  day  and  location 

a  specific  volume 

a.s  scattering  albedo  for  the  atmosphere  (if  clouds  are  present 

they  determine  the  scattering  aibedo) 
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(Xg  albedo  of  ground  surface  as  a  function  of  hour  angle  and 

surface  characteristics  (Wetzel  1978) 

V  horizontal  gradient  operator  (=  i3/3x  +  j8/9y) 

7  potential  temperature  lapse  rate  above  the  PBL 

A0  jump  discontinuity  at  the  top  of  the  boundary  layer 

Au  vector  difference  between  the  wind  just  above  the  boundary 

layer  (V0)  and  the  wind  in  the  boundary  layer 

(Au  =  Atii  -f-  Atij) 

0-  potential  temperature  at  the  top  of  the  PBL 

0O  potential  temperature  at  the  ground 

0i  potential  temperature  of  the  free  atmosphere  above  the  PBL 

0P  potential  temperature  above  jump  discontinuity 

0®  potential  temperature  at  the  top  of  surface  layer 

0„  potential  temperature  at  the  top  of  the  PBL 

0*  — SH /ii* 

0  mean  potential  temperature 

T)  boundary-layer  coordinate  in  vertical  direction 

V  boundary-layer  coordinate  vertical  velocity 

K  thermal  conductivity  of  soil 

X  latent  heat  of  evaporation 

U  kum/fL 

Mi  h/L 

£  angle  between  the  surface  stress  and  the  geostrophic  wind 

above  the  boundary  layer 

P  density  of  air 

Pvj  density  of  water  (1  gm  cm-3) 

o  Stephan-Boltzmann  constant 

o >,  height  of  the  boundary  layer  in  a-coordinates 

tp  period  of  cycle 

t,  t  mean  total  transmission  functions  for  effective  absorber  u*  at 

temperature  T 

Txy  surface  stress  vector 

<t>  geopotential 

(P  frequency  of  the  variance  of  soil  surface  temperature 

uJbV  Brunt-Vaisala  frequency  =  (^) 
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APPENDIX  B 


Derivation  and  Energy  Conservation  of  the  Eta-Coordinate  System 


a.  Derivation  of  the  equations  in  ^-coordinates 

The  most  common  form  of  terrain  following  coordinate  is  the  cr-system  where 


a 


(P 


-  Pt) 

7 


(B.l) 


where  7  —  (ps  —  pt),  and  p®  and  pt  are  the  pressures  at  the  surface  and  the  top  of 
the  model  domain,  respectively.  The  equation  set  in  this  coordinate  system  may  be 
written  (Haltiner  and  Williams  1980) 


^  +  V0  +  oaVx  -f  fk  X  V  -  F 

(B.2) 

a 

I 

II 

■*|b 
(0  1(0 

(B.3) 

If  +  v  'v  +  *ti  -  o 

(B.4) 

-  a"  -  « 

(B.5) 

where 

and 


IU  =  ~  =  TO  4-  0(97-  /dt  -f  V  ■  V 7T ) 
at 

V  =  i9/9x  +  j9/9y 


In  order  to  preserve  conservative  properties  in  the  finite  difference 
equations,  we  write  the  total  derivatives  in  flux  form,  which  can  be  obtained  with 
the  aid  of  the  continuity  equation  (B.4).  Equations  (B.2)  and  (B.5)  can  then  be 
written 


drV  ,  durV  ,  dvrV  ,  SxVcr 
dt  dx  dy  do 


—  —  r(xoVr  —  / k  X  irV  +  tF 

(.B.6) 
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(B.7) 


drT  ,  BurT  ,  BvkT  ,  dr  Ter  _  r<x.  . 

3f  "r  3x  ^  By  ■r  3<r  cPw  ^  cP 


The  transformation  to  T7-coordinates  is  carried  out  formally  by  using  the 
definition  of  V 


H 


Ok 

1  ~~Oh 


so  that 

a  =  rjH+a^ 


cr  <  oh 
cr  >  a* 


(B.8) 


(B.9) 


and  for  some  dependent  variable  A 


nyi  =  i  <y 

IScrJcr  H  377 


f<y]  _  fa^41  I  BA  (3a) 

l3sV  =  ^Bstrj  H  377  l3s>r? 


where  s  —  x,  y,  or  f.  Then  the  H-coordinate  set  equivalent  to  (B.6),  (B.3),  (B.4),  and 
(B.7)  is 


BrHV  ,  BurHV  BvrHV  ,  BrWHV 

3 1  3x  By  +  BV~ 

-  rHV<t>  -  rHcxVar  -  f k  x  tHV  4-  rH F 


-rHa 


rHV  +  =  0 

a7! 

BurHT  BvrHT  BrTHy 

Bx  By  377 


TrHa, ,  , 

<*>  +  -f— 


where  now 


uj  =  =  rHf]  -+-  3(airj/3f  -f  V  •  Vctt 

at 


(B.10) 

(B.ll) 

(B.12) 

(B.13) 

(B.14) 


Note  that  a  is  now  a  dependent  variable  which  is  a  function  of  x,  y,  V  and  (,  and  is 
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given  by  (B.9).  As  can  be  seen  by  comparing  (B.10)-(B.13)  with  (B.6),  (B.3J,  (B.4),  and 
(B.7),  the  transformation  leads  to  the  prognostic  variables  being  weighted  by  irH  in 
the  V  system  instead  of  the  ir-weighting  present  in  the  a  system.  Note  also  that 
terms  that  had  been  a  times  derivatives  of  ir  become  derivatives  of  the  quantity 
fair)  in  the  T)  system. 

b.  Energy  conservation  in  the  H-system 

In  order  to  determine  the  correct  finite  differencing  form  that  will  perserve 
energy  conservation  in  the  finite  difference  equations,  the  energy  conservation 
constraints  of  the  continuous  equations  must  first  be  derived.  The  analysis 
presented  here  closely  follows  that  carried  out  by  Haltiner  and  Williams  (1980)  for 
the  a-coordinate  equations. 

We  begin  with  the  ^-coordinate  momentum  equation  [equivalent  to  (B .2)1 

^  +  V*  +  aVair  +  /k  >  V  =  F  (B.15) 

at 

This  is  dotted  with  tHV,  and  rewritten  in  the  flux  form  with  the  aid  of  the 
continuity  equation  (B.12),  to  yield 

+  V  •  (ItHVV2)  +  JL(i*H77t'2)  = 

-irf/V  (tHV*  -f  a Vair)  +  r H\  F  (B.16) 

The  first  term  on  the  right  hand  side  represents  the  kinetic  energy  production  by 
the  pressure  force.  We  expand  this  term  as  follows 

-irHV  •  (V*  4-  aVair)  =  -V-  ( irWVtf)  +  *V-  firWV)  -  air HV  •  Vair 

The  continuity  equation  (B.12)  allows  this  to  be  rewritten 

=  -V  •  (ir H\<t>)  -  -  air HV  ■  Vair 

=  —  V  .  +  *Hr,%  ~  aTHV  '  Vcrr 
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Use  of  the  hydrostatic  equation  (B.ll)  leads  to 


=  -V  •  (tHV^)  -  +  *HTH-Htcl)  -  a tHV  •  Vctt 

dt  0T7 

Adding  and  subtracting  t Ha(daK  /dt),  and  rearranging  yields 

-irHV  •  (V*  +  aVar)  =  -V  ■  (irHV*)  -  +  rHa^ 


-  Hi rct(^  +  V  ■  Vcttt  +  irtfr?) 


Then  using  the  definition  of  u),  (B.14),  and  expanding  the  time  derivatives  of  the 
products  leads  to 


=  -V  •  (ir H\<t>)  - 


30ir  HT] 


~  *’«<* if 


—  tHclu) 

Now,  (B.ll)  can  be  rewritten  as  d{j>o) /d*)  —  —  B(ircra  —  <t>),  and  this  along  with  use 
of  (B.ll)  directly  leads  to 

_  n  ,.uv^  3 0tHV  971  _3tf  9ct 

-  -V  •  (*HV*)  -  -K—  —  -  T/yaw 

Noting  that  t  is  not  a  function  of  Tj  and  that  3(3cr  /dt)  /dT)  =  dH /dt  leads  to  the 
result  that 


—  irHV  •  (V*  4-  aVcnrt  =  —  V-  +  £cx|y) 

-  -  THaW 


(B.l  7) 


We  form  the  total  energy  equation  by  using  (B.16)  with  the  substitution 
given  by  (B.l 7),  and  adding  it  to  the  thermodynamic  equation  (B.13)  which  can  be 
written  in  the  form 


|f(THcPT)  4-  V.  ixHVcpT)  4  ^vHcpTHfi)  =  rH{a w  +  Q) 
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(B.18) 


to  yield 


<HV2  +  tHCpT)  +  V-  V(iir HV2  +  -kHc^T  +  tH*) 

+  UiTVlHri  +  Tc-7Hil  +  ’♦H  +  M*°W  +  |kf ) 


=  *H«3  +  V  ■  F) 


(B.19) 


If  (B.19)  is  integrated  from  T)  =  —  1  (where  cr  =  0)  to  r?  =  1  (where  a  =  1),  we 
obtain 


4-  J  irHfiv'2  +  CpTjdr? 


-L 


rH(Q  4-  V  •  F)dr? 


CpT  4-  0jd57 
(B.20) 


where  we  have  used  Hfj  =  0  at  r?  =  —1  and  1,  3cr/at  —  0  at  r?  =  —  1  and  1,  and 
$(3ir /9 1)  —  d(6sps)  /dt.  This  is  precisely  the  same  energy  relation  derived  by 
Haltiner  and  Williams  [1980,  their  Eq.  (7-42)]  for  the  a-coordinate  equations,  if  we 
note  that  HdT)  =  da.  It  is  desirable  for  this  energy  constraint  on  the  continuous 
equations  to  also  hold  for  the  finite  difference  equations.  As  will  be  seen,  this  will 
determine  the  appropriate  form  for  the  finite  difference  equation  set,  and  for  the 
method  of  vertical  finite  differencing. 


c.  The  energy -conserving  finite  dif f erence  equations 
The  finite  difference  form  of  (B.10)  is 

+  A(u**HVk)  4  ^(vkr HVJ  + 

Hlk-t-l  /2^\  +  1  /2)  —  (H7))*  -  I /2Vk_j /2)J  =  (B.21  ) 

—  4-  akV(7vr]  —  / k  v  rHVk  rHFk 

where  H  takes  on  the  appropriate  value  as  given  by  (B.8)  depending  on  whether  the 
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layer  fc  is  above  or  below  <jh.  The  indices  (fc  — 1/2)  and  (k+1/2)  refer  to  the  layer 
interfaces  bounding  layer  fc  above  and  below,  respectively.  The  quantity  V  is  the 
interpolation  of  V  from  adjacent  layers  to  the  interface.  Haltiner  and  Williams 
(1980)  show  that  the  advective  terms  will  conserve  both  V  and  (V  V)/2  if 


v*+1/2  =  ^(Vk+1  +  vk) 


(B.22) 


We  wish  to  insure  energy  conservation  of  the  rhs  as  well.  The  Coriolis  terms  do 
not  contribute,  and  we  will  ignore  the  friction  term  and  concentrate  on  the  pressure 
gradient  term.  As  with  the  continuous  equations,  we  find  the  rate  of  working  by 
the  pressure  gradient  terms  by  taking  Vk  •  of  the  first  term  on  the  rhs  of  (B.21)  to 
obtain 


-TtfV*-  [V*K  4-  =  -V-  (* HVk*k)  +  0kV-  (r HVk) 


—  CXk-K HVk  ■  Vct*  7T 


(B.23) 


The  finite  difference  form  of  the  continuity  equation  is 

^  +  V-  tH\k  +  ^[(Hi7)k  +  l/i  -  \Hrj\-  \/<j  =  0  (B.24) 

so,  (B.23)  can  be  rewritten 

—  THVy  [V0k  -f  a^Vcr,.7r]  =  —V  -  (TrH\k<t>k) 

Adding  and  subtracting  ^['H^)k+1/20k+1/2  -  (H77)k_1/20k_1/2]  yields 

=  -V  -  !irHVk0fc)  -  0*^  -  ^[(H4  +  ./^  +  ./i  -  (H»7)fc-i/^*-i/:] 

+  T^[(ff77;*.  +  i/2(0k  +  i/i  ~~  (Hr?'k_i/2(0k_)/;.  —  0k 


—  tHV,  ■  Vakr 
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Now,  adding  and  subtracting  Hzak(dakir  /dt)  leads  to 


—  *H\k  ■  [V*k  4-  a*VCTk*]  =  -V-  (t HVk<(>k) 

~  1/20*1- 1/2  —  (f/77)k_)/20k_i/2J 

-  {4>kH  -  H*akc7k)^  -  ^x|  +  H*2ak^ 

(B.2S) 

-  +  V*.  Vakrj 

+  ^[(Hl7)k+l/2(^k+i/2  —  0k)  —  (H??)k_1/;.!0k_i/2  —  0k)J 


Now,  it  is  easy  to  verify  that 


3cr*  =  OkdH 

dt  h  dt 


(B.26) 


So,  (B.2S)  can  be  written 

*HVk  •  [V<>k  +  akVcr*ir]  =  -V-  (t HVk<>k) 

~~  £^[{Wi7)k+i/20k+i/2  —  (W77)k_1/20jt_)/2J 

-  [<t>kH  -  HTrakCTt]|5  -  jj[*kH  - 

(B.27) 

-  akTTHp^  4  Vk  •  VcrkT 

+  ^[(H’7)lc  +  1/2(^k+1/2  —  0k)  —  (HT7)k_i/2(0k-i/2  —  0k)J 


But,  using  the  finite  difference  form  of  d{^cr)/dT]  =  —Hiraa.  —  0),  and  noting  that 
d* /dt  and  dH  /dt  are  not  functions  of  7?,  we  can  write 
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—  irHVk  [V0k  altVCTkir]  =  -V- 


7)k  +  i/2^k  +  l/i  —  (W?7)k-l/20k-l/2j 


-  ^^»c  +  l/2^k  +  l/2^  —  ^k-l/2^>c-l/2^J 

1  h  _^k  +  l/2  3H  a  .^k-l/23H] 

-  H—3F  “ 

-  a^w[^  +  Vk  •  Vakir] 


+  ^^^H?7)k+1/2(^k  +  i/i  ~  0k)  —  (HJ7)t_, 


(B.28) 


/2\Vk-l/2 


*k>] 


Using  (B.26)  again  leads  to  our  final  result 


*HVk-  [V0*  -  akVCTk?r]  >=  -V-  (irHVk0k) 


lCTk  +  l/20k  +  l/2^j-  —  <7k-l/2^k-l/2 


- 

Ar?[ 

1  II  k  +  1/2  -  _9<5-k-i/2 

~*t  0ic-'/^~ dj— _ 

•] 


a.kirH 


dakr 


—  +  Vk  •  V  a  k  ic 


(B.29) 


+  ^^[(H^)k+l/2(^k+l/2  —  0k)  —  (HV)k- l /?{$*. -1/2  —  0*)J 


On  term  by  term  comparison  with  (B.17),  it  is  clear  that  (B.29)  will  match  the 
continuous  equations,  and  conserve  energy  in  the  same  way  on  summation  over  all 
layers,  if  we  specify  that  w  is  defined  by 

(B.301 

(HH)k- i/2(ik-i/i  -  0kf| 


rH  a 


:ujk  =  akir  H 


3£V 

3f 


-t  V*  ■  V<7k* 


+  ,  /2(0k+  l ,1  —  0k)  — 
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This  definition  of  U  is  precisely  what  is  needed  to  specify  the  finite  difference 
form  of  the  thermodynamic  and  hydrostatic  equations,  which  we  will  derive  next. 

In  order  to  show  the  constraints  placed  on  the  finite  difference  form  of  the 
thermodynamic  equation,  we  need  to  first  show  two  alternative  finite  difference 
forms  of  the  total  derivitive.  For  a  variable  A  at  level  fc,  we  can  write  the  total 
derivative  in  “flux  form’’  as 


+  V*  •  VrHAk 


(B.31 ) 


T  PJk+l/Z-^k  +  l/?  l/2.4k-l/2J 


But  use  of  the  continuity  equation  (B.24)  allows  this  to  be  rewritten  in  an 
“advective  form”  consistent  with  the  flux  form  as 


rH  lx  +  Vk  VAk  + 


(B.32) 


:[(F/t7)k  +  1/2(i4»;  +  1/2  —  /Ik)  —  (Ff77)k-i/2i-4k- 1//  — 


For  frictionless,  adiabatic  motion,  the  potential  temperature  is  conserved,  so 
dQ  / dt  =  0.  The  finite  difference  version  of  this  in  flux  form  is 

+  V  •  THVk6k  +  ^[(WW./A+i/:  -  (W77)k-1/20k-i/J  =  0 


(B.33) 


In  order  to  have  0  and  9 2  conserved  (Haltiner  and  Williams,  1980),  we  define 


0fc  +  l/2  —  ^ 


(B.34) 


We  take  the  potential  temperature  to  be  defined  by 


0k  =  Tk/Pk 


(B.35) 


Pk  —  j^(Pk-l/J  +  Pk  +  l/2)  / lOOoj* 


(B.36) 
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The  conservation  of  &  can  be  written  in  advective  form,  (B.32),  as 


th[^  +  V,  VGk  + 

^^[(H?7)k+i,/2(0j.+i/2  —  8k',  +  (Hfj)k_,/2(6k  —  9k_J/2)J  =  0 

Substituting  the  definition  of  0k,  (B.35),  this  may  be  rewritten 

*»[!  + v*  -  vlr‘  +  +  v*  ■  (B-37) 

^^(Ht7)k  +  l/2(Pk0k.  +  I/2  —  Til)  +  {HV)k.-l/2,'.Tk  —  p  k0k_i/2)J  =  0 

where  we  have  made  use  of  the  fact  that  Pk  can  be  considered  a  function  of  a  and 
it.  We  now  add  Tk  times  the  continuity  equation,  (B.24),  so  that  the  first  term  can 
be  rewritten  in  flux  form,  add  the  finite  difference  form  of  ir9T/3^  to  both  sides, 
and  multiply  the  whole  equation  by  cp  to  obtain 

~ CpitHTj.)  +  V  •  (cPirHTkVk)  +  C^[(Hr?)k+1/ ,Tk+!/2  +  (Hr7)k_,/2tk_1/2)] 

-  +  v*  •  lB'38> 

+  ■^y[(Hn)k  +  l/2(’Pk-kl/2  —  Pfc0»:  +  l/2)  +  (HT7)k_1/2(P  k0k_i/2  —  Tk_1/i)j 

The  lhs  of  (B.38)  is  the  finite  difference  form  of  the  lhs  of  (B.18),  and  therefore, 
the  rhs  should  be  equal  to  irHaw  (since  we  are  assuming  Q  =  0  here).  Comparison 
with  (B.37)  shows  that  the  first  terms  will  be  equal  if 


__  CpT k  dP y 
Pk  d(7kr 


(B.39) 


Equating  the  other  terms  leads  to  the  following  relations 


Cp(Tk  +  i  —  P  k0k+i)  —  4>k  —  0k  +  i 
CpiP  k0k_i  Tk_i)  =  4>k~\  4>k 


(B.40) 
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When  (B.35)  is  used  in  these,  they  can  be  rewritten 


(CpT'k  + 1/2  +  ^fc+1/2)  —  (CpTk  +  4k)  =*  P  kCp(6k.  +  !/2  —  0fc) 
(CpTfc  +  4k)  —  (CpT*._1/;  -f  ^k-1/2)  “  P  kCp(@k  —  ®k+J/z) 


(B.41) 


We  can  replace  fc  by  fc  +  l  in  the  second  of  (B.41),  add  it  to  the  first,  and  again  use 
(B.35)  to  obtain 

4k+ 1  —  4k  =  — Cp(Pk  + 1  —  Pk)&k+ 1/2  (B.42) 


This  represents  a  finite  difference  form  of  the  hydrostatic  equation  that  is 
consistent  with  the  other  equations,  and  provides  a  means  of  calculating  the 
geopotentials  at  all  layers  once  the  geopotential  is  known  at  the  lowest  layer  (where 
k  =  kbm).  Haltiner  and  Williams  (1980)  show  that  it  is  possible  to  derive  an 
integral  constraint  from  the  energy  conservation  which  yields  the  geopotential  of 
the  lowest  layer  in  terms  of  the  geopotential  differences  of  all  the  other  layers.  It 
is  pointed  out,  however,  that  this  accumulates  the  errors  of  the  layer  calculations 
and  leads  to  large  errors  in  the  geopotential  of  the  lowest  layer.  Experimentation 
with  the  model  verified  this  result.  We  choose  a  simpler  method  of  obtaining  the 
lowest  layer  geopotential  which,  though  not  strictly  consistent  with  energy 
conservation,  yields  accurate  values.  Since  the  boundary  layer  parameterization 
provides  an  “surface”  temperature,  we  use  this  temperature  to  form  a  “half-layer" 
average  potential  temperature 


^kbm  +  \  /*  2^k*>m  ® s  f  c) 


Then  an  equation  of  the  form  of  (B.42)  can  be  used  to  find  4kbm  in  terms  of  the 
geopotential  at  the  surface. 

The  thermodynamic  equation  can  now  be  written  in  the  form  used  in  the 
model  by  noting  that  (B.38)  reduces  to 


dr  HTk 
dt 


+  V  •  *HVkT* 


r_ 

&rj\ 


[(HWk  +  i 


/2  P  k&k  +  \, 


,'.P'h)k-  \  nP k 0 


7k-l 


•J 


•*HaJ  3 
cp  dt 
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(B.43) 


(where  we  have  added  the  diabatic  heating  term).  Thus,  the  equations  (B.21),  (B.42), 
and  (B.43)  form  a  consistent  finite  difference  set  for  momentum,  geopotential,  and 
temperature  which  conserve  total  energy  in  the  same  manner  as  the  continuous 
equations. 

d.  The  complete  finite  dif ference  equation  set 

We  restate  here  the  finite  difference  equations  derived  above  and  add  those 
not  yet  discussed  to  form  the  complete  set  used  in  the  model.  The  only  advective 
quantity  in  the  model  whose  finite  difference  form  is  not  discussed  above  is  the  one 
that  governs  specific  humidity,  q.  In  the  absence  of  condensation  or  evaporation, 
however,  specific  humidity  represents  a  conserved  quantity  that  satisfies  dq/dt  = 
0.  We  write  this  total  derivative  in  the  flux  form  given  by  (B.31). 

The  complete  set  of  finite  difference  equations  for  momentum,  4>,  T,  and  q  is 

then 


^[(W77)k+i/2V*-n/2  :  -  (HJ7)k-i/2Vfc_1/2)]  -  (B.21) 
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where 


(B.35) 


Qk  =  Tk/Pk 

and 

Pk  =  [|(p*-„2  Hr  Pk+.^/lOOoj*  (B.36) 

Note  that  an  eddy  diffusion  term,  Fk,  is  added  to  the  rhs  of  (B.43')  and  (B.44)  as 
described  in  section  2  to  help  control  noise.  The  finite  difference  form  of  (2.7), 
which  calculates  the  rate  of  change  of  r  is  written 

Kbtn  -. 

If  -  -  (B-45) 

where  H  takes  on  the  correct  values  above  and  below  the  boundary  layer  top  as 
given  by  (B.8).  Similarly,  the  finite  difference  form  of  (2.8)  allows  calculation  of 
(H?})  at  each  interface  as 

(HT})k+i/2  —  —(r]k+i/2  +  l)  +  ^IfJ 

(B.46) 

-  life"™*'  + 

which  can  be  applied  after  (B.45)  has  been  solved  to  provide  3  r /dt  and  after  the 
boundary  layer  parameterization  has  provided  dh/dt  which  can  be  converted  to 
dah/dt.  The  above  finite  difference  equations  are  applied  on  the  staggered  grid 
using  the  averaging  and  differencing  schemes  presented  by  Anthes  and  Warner 
(1978). 
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APPENDIX  C 


Creation  of  Terrain  Fields 

a.  Creation  of  unmodified  and  envelope  terrain 

This  appendix  describes  the  procedure  followed  in  developing  the  terrain 
fields  used  for  the  simulations  discussed  in  section  4.  A  terrain  field  is  simply  an 

array  that  has  an  elevation  value  for  each  grid  point,  but  unless  care  is  taken  in 

specifying  these  elevation  values,  the  terrain  field  can  generate  spurious  noise  in 
the  model — either  because  of  an  amplifying  resonance-type  response  at  2£z 
wavelengths  or  because  of  inconsistencies  between  the  CGM  and  FGM  terrain  fields 
in  the  region  of  overlap  between  the  two  meshes. 

As  an  inital  step  in  creating  the  terrain  field  used  here,  elevations  were  read 
off  U.S.  Geologic  Survey  topographic  maps  at  the  latitude/longitude  locations  of  the 
FGM  and  CGM  grid  points.  In  addition,  with  each  grid  point  treated  as  the  center 
of  a  20  km  (FGM)  or  60  km  (CGM)  grid  box,  the  elevation  of  eight  additional  points 

that  were  within  the  grid  box  were  read  off  the  maps — so  that  the  model  grid  point 

was  the  center  point  of  a  3X3  array  of  equally  spaced  points  in  the  grid  box.  For 
each  grid  box,  the  nine  points  were  used  to  create  and  average  elevation  and  a 
standard  deviation  of  elevations  within  that  grid  box.  The  average  plus  the 
standard  deviation  for  each  grid  point  produces  a  terrain  height  field  referred  to  as 
a  partial  envelope  orography  (Wallace  et  al.  1983).  For  a  full  envelope,  the 
standard  deviation  would  be  multiplied  by  2.0  (referred  to  as  the  envelope 
parameter).  An  envelope  orography  is  desirable  because  it  is  better  able  to  produce 
the  blocking  effect  of  mountainous  terrain,  which  is  underestimated  by  a  simple 
average  over  the  grid  box. 

b.  Smoothing  and  modification  to  ensure  CGM/FGM  compatibility 

While  the  use  of  the  envelope  technique  accomplishes  some  smoothing  of  the 
terrain  fields,  it  is  necessary  to  cary  out  some  additional  filtering  to  remove  strong 
harmonics  at  2Ax  wavelengths.  This  was  accomplished  by  applying  a  two- 
dimensional  Shapiro  (1970)  filter  on  the  terrain  field.  In  order  to  avoid  the  need  for 
a  modified  filter  near  the  boundaries,  the  envelope  terrain  field  was  initially  created 


with  an  extra  row  of  grid  points  on  all  sides  of  the  desired  final  domain  size.  After 
the  filter  had  been  applied,  these  extra  points  were  removed  from  the  array  to  yield 
the  smoothed,  partial  envelope  CGM  and  FGM  domains  at  the  size  required  for  the 
model. 

The  last  step  in  the  processing  involved  a  modification  to  the  FGM  grid  point 
elevations  near  the  FGM/CGM  interface.  Zhang  et  al.  (1986)  have  pointed  out  the 
strong  need  for  consistency  in  the  CGM/FGM  overlap  region  (see  Fig.  2)  to  help 
reduce  noise  generation,  especially  since  the  two-way  interactive  boundary  condition 
results  in  an  overspecification  of  the  variables  here.  In  the  two-way  interactive 
scheme,  the  FGM  solution  is  used  at  CGM  boundary  points  after  the  application  of 
a  nine-point  average  of  the  FGM  variable  centered  on  the  CGM  grid  point.  In  order 
for  these  averaged  FGM  quantities  to  be  consistent  with  the  appropriate  CGM 
quantity,  the  average  FGM  elevation  using  the  same  averaging  operator  must  equal 
the  CGM  elevation  at  the  coincident  point.  Zhang  et  al.  (1986)  describe  a  method 
for  modifying  the  terrain  field  so  that  this  condition  is  satisfied,  and  their  method 
was  employed  on  the  terrain  field  used  here.  The  resulting  terain  fields  after  all 
this  processing  are  those  shown  in  Figs.  15  and  18. 
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